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FOREWORD 

The Mark XV Error Propagation Program was developed under Contract 
NAS1-9307 for the National Aeronautics and Space Administration, Langley 
Research Center by Philco-Ford Corporation. Two tasks were undertaken 
in this development. They were: 

1. Conversion of the Mark II Error Propagation Program to 
a form acceptable to the CDC 6600 computer, and 

2. Addition of capabilities for start-up and search for 
suitable interplanetary trajectories. 

The work was performed in the period June 30, 1969 through December 31, 
1969. The Philco-Ford personnel responsible for the program development 
and writing of this manual ware: 

W. S. Bjorkman Senior Engineering Specialist 

M. J. Brooks Senior Programmer 

The effort was managed by R. C. Jensen. 

This manual contains a description of the Mark IV Error Propagation 
Program and detailed instructions for its utilization. For additional 
details about the implemented theory, refer to "User's Manual for the 
Mark II Error Propagation Program", WDL-TR2758, and to "Subroutine 
Descriptions and Listings for the Mark II Error Propagation Program", 
WDL-TR2757, Volume I. 

It is Philco-Ford' s hope and expectation that the Mark IV program will 
prove helpful in- understanding the interrelationships of the various 
parameters which affect mission success and will, therefore, become a 
useful tool in the planning of future space missions. 
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INTRODUCTION 

Space mission planning requires an understanding of the interrelationships 
among the numerous parameters affecting mission success. In most instances, 
no simple equation defining these interrelationships exists. The only 
practical way of studying the effects of many of these parameters is to 
simulate the mission on a digital computer, using the parameters together 
with laws of physics and principles of statistics in the simultation. 

The Mark IV Error Propagation Program is a computer program which enables 
the study of parameter interrelationships in the realms of lunar or 
interplanetary trajectories, navigation and guidance. It is a fourth- 
generation program which was developed for LRC by Philco-Ford. Three 
earlier error propagation programs were developed by Philco-Ford under 
NASA contracts for GSFC. The Mark IV program is a version of the GSFC 
Mark II Error Propagation Program converted for operation in Fortran IV 
on the CDC 6600 computer and modified by the addition of start-up and 
search capabilities. The Mark II program, delivered in December 1965, 
was written partly in Fortran IV and partly in machine language (MAP) 
for the IBM 7094 computer. It lacked the start-up and generalized search 
capabilities, presuming the initial trajectory conditions to be obtained 
from some other program. A restricted search capability for integrated 
trajectories was available in the Mark II program, however. The Mark IV 
program retains all of the capability of the Mark II program, is gen- 
erally more efficient, and has been simplified in its programming. 

This manual is divided into two sections. The first contains a brief 
description of the capabilities of the Mark IV program and includes a 
summary of the theory, assumptions and equations which have been 
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implemented. The second section contains instructions for using the 
program along with input and output examples. Tables of user-relevant 
information are provided for easy reference once a program familiarity 
has been developed. The user is referred to the user’s manual (Ref. 1) 
and subroutine descriptions (Ref. 2) for the Mark IX program for addi- 
tional details not found in this manual. 
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SECTION I 

PROGRAM CAPABILITIES 

The capabilities of the Mark IV Error Propagation Program fall into two 
major categories which are: 

1. trajectory generation, and 

2. error propagation. 

The error propagation capability requires a nominal trajectory which must 
be generated and stored on tape or disc sometime prior to propagation of 
errors. We will, therefore, begin this section with a description of the 
program's capabilities for trajectory generation and conclude with error 
propagation capabilities. 


1.1 TRAJECTORY GENERATION 


The sections of the program concerned with trajectory generation are 
called START-UP, SEARCH, CONW and PINT. START-UP provides approximate 
initial conditions for interplanetary transfer trajectories using the 
matched-conic assumption. SEARCH is a "generalized" search routine which 
iterates by the steepest descent method to differentially correct initial 
conditions in order to satisfy imposed constraints at the end of the 
trajectory. The trajectory model used by SEARCH may be chosen to be 
either patched-conic or precision-integrated. CONW takes specified 
initial conditions for a patched-conic trajectory and writes interpola- 
tion coefficients for that trajectory on tape or disc for later use in 
error propagation. PINT'S primary function is to write integrated 
trajectory interpolation coefficients on tape just as CONW does for 
patched-conic trajectories. But, PINT can also perform a restricted 
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search on initial conditions to satisfy end constraints as SEARCH does, 
using an integrated trajectory model. Another PINT capability is to 
compute and integrate variational equations for equation of motion 
parameter sensitivities or the state transition matrix. CONW and PINT 
were sections of the Mark II Error Propagation Program and will, therefore, 
be discussed only briefly. 

1.1.1 START-UP 

This capability requires the user to specify the launch and target planets 
plus the departure and arrival dates for interplanetary transfer. The 
program then interrogates the planetary ephemeris for positions and 
velocities of the launch body at departure date and target body at arrival 
date. The heliocentric conic section which joins the two (massless) 
planets in the specified flight time is then determined. The vector 
difference between the initial velocity on the heliocentric conic and 
the launch body velocity at departure date is taken to be the hyperbolic . 
excess velocity at launch. 

The next assumption used in the determination of approximate initial 
conditions is that the trajectory originates from a circular parking 
orbit about the launch body. The user must specify four parameters which 
describe this parking orbit, namely: insertion latitude, insertion longi- 
tude, insertion velocity azimuth and orbit altitude. If the launch body 
is the Earth, the START-UP assumes the parking orbit parameters to be 
geographic (Earth-fixed) and uses them to find the time of day at which 
the hyperbolic excess velocity vector lies in the plane of the parking 
orbit. If the latitude of the hyperbolic excess velocity vector is 
greater than the inclination of the parking orbit, the program finds the 
time of day which minimizes the latitude of the hyperbolic excess velocity 
vector relative to the parking orbit. The time of day thus found is 
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interpreted to be the time of insertion into parking orbit (new departure 
date) and is always later than the input departure date. If two solutions 
for time of day exist, the one nearer the input departure date is selected 
by the program. If the launch body is not the Earth, the parking orbit 
parameters are assumed to represent an “inertial" parking orbit which does 
not rotate with the launch body. The program next proceeds to calculate 
the time in parking orbit which results in the minimum injection velocity 
requirement and yet attains the desired hyperbolic excess velocity. The 
injection maneuver is a velocity impulse. The quantities output from 
(i.e., provided by) START-UP for use by SEARCH, CONW or PINT are: 

Date (and time, UHE) of park orbit insertion 

Date (and time, UMT) of injection onto the transfer hyperbola 

Launch parameters - park time, azimuth of the velocity impulse 

relative to the park orbit track at injection, 
elevation of the velocity impulse at injection, 
and velocity impulse magnitude 

Cartesian state - three components of position and three of 

velocity at injection referred to the mean 
equator and equinox of 1950.0 coordinate system 

Spherical state - radius, declination, right ascension, speed, 

flight path angle and azimuth at injection 
referred to'EESO. 

* * 

The conditions provided by START-UP are precise enough to initiate a 
patched-conic or integrated trajectory to the vicinity of the target body.- 
An option for “matching'' conics is provided if a better solution is 
desired. Under this option, the program utilizes the previously-computed 
departure trajectory to calculate the point and time of “patch 1 ' to the 
sun. The hyperbolic excess velocity and specified desired miss vector 
at arrival are used to compute the arrival hyperbola and thus the point 
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and time of "patch" to the target body. The heliocentric trajectory is 
then re-computed as the conic section which transfers between the sun- 
patch and target -patch points in the adjusted transfer time. Of .course, 
the launch and target bodies move during the departure and arrival 
phases. The new heliocentric conic results in new hyperbolic excess 
velocities at departure and arrival which may then be used to initiate 
another iteration. Otherwise, the hyperbolic excess velocity at de- 
parture is simply used to provide an improved set of launch control 
parameters and injection conditions according to the parking orbit 
assumptions. A summary of the mathematical formulation of START-UP 
may be found in Appendix A. 

1.1.2 SEARCH 


This program performs a systematic search to satisfy a specified set of 
end constraints by differentially correcting a set of initial control 
parameters. The control parameters (from one to six in number) may be 
selected from any one of the following sets. 


Set 1 Set 2 


Set 3 


Cartesian 

X 

Y 

Z 

vx 

VY 

vz 


Spherical 

R 

LAT 

LON 

V 

Y 
AZ 


Launch Parameters 

DIAZ (insertion azimuth) 

DTL (time of insertion) 

PRKT (park time) 

DVAZ (injection impulse azimuth) 
DVEL (injection impluse elevation) 
DELV (injection impulse magnitude) 


The cartesian and spherical controls may be referred to the Earth's 
equator and equinox or to the ecliptic and equinox (both mean of either 
1950.0 epoch or date). Their origin may be a time other than when 
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initial conditions are specified, making it possible to target from a 
midcourse maneuver time. Choice of the launch parameter set as control 
variables requires that the parking orbit and insertion date be speci- 
fied and that the search take place at injection from parking orbit into 
the transfer orbit. 

The constraints may be selected from: 


1. 

- 6. 

cartesian end state components 

7. 

-12. 

spherical end state components 

13. 

B.T. 

, asymptotic miss vector component 

14. 

B.R. 

, asymptotic miss vector component 

15. 

V oo > 

hyperbolic excess speed 

16. 

V 

radius at periapsis 

17. 

i. 

inclination 

18. 

n 

, longitude of the ascending node 

19. 

V 

argument of perifocus 

20. 

V 

time of flight 


The number of constraints cannot exceed the number of controls and is 
always six or less. The frame-relative constraints may be referred to: 

1. Earth's mean equator and equinox of 1950.0, 

2. target orbital coordinates (1. along the radius vector from 
the target's central body, 2. normal to the other axes in 

the right-handed sense, and 3. along the target's orbit normal), 

3. targetographic coordinates, or 

4. mean ecliptic and equinox of date. 
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The trajectory model (or plant) by which constraint errors are computed 
from controls may be either patched-conic or precision-integrated to 
include perturbations. In the integrated formulation, Encke's Method 
is used for trajectory calculation and Adams' Fourth-Order Method is 
used for numerical integration. Sensitivities of constraints to controls 
are computed by the secant or difference method. That is, finite incre- 
ments are added to each control in turn and the trajectory and constraints 
are re-calculated. The sensitivities are then taken to be the constraint 
differences divided by the control increment. This technique is costly, 
but reliable and conceptually simple. The Method of Steepest Descent is 
used to predict control increments which will reduce the constraint 
errors. 



Figure 1-1 Search Logic 

Figure 1-1 is a block diagram of the computational logic of the 
SEARCH program. The symbol Y represents the computed constraint vector 
and T ^ the desired constraint vector. 
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A scanning option has been provided in SEARCH whereby any control variable 
may be automatically incremented. This option is useful for assessing the 
nature of local constraint behavior due to control variation. 

Appendix B contains a brief description of the mathematical formulation 
of SEARCH. A description of the acceleration equations used in the 
integrated trajectory model is found in Appendix B of Reference 1. 

Appendix C of this manual describes the method of numerical integration 
used here. 

1.1.3 CONW 

This section of the Mark IV program accepts initial trajectory conditions 
and computes a patched-conic trajectory. Interpolation coefficients for 
the trajectory are written on tape or disc for later use by ERP, the 
error propagation section. In addition to trajectory interpolation 
coefficients, CONW computes the matrix of sensitivity of end conditions 
(B.T, B.R, time to periapsis, v^ , r^, i) to end cartesian state for 
later use in guidance and prediction calculations. The advantage of 
using a patched-conic nominal trajectory lies in the speed and simpli- 
city of its generation - either type trajectory gives rise to similar 
answers in error propagation. Initial conditions for CONW may be ob- 
tained from START-UP, SEARCH or other sources. CONW accepts cartesian 
or spherical components of state or orbital elements referred to equa- 
torial, ecliptic or body-fixed coordinates of 1950.0 or date. 

1.1.4 PINT 

This section of the Mark IV program accepts initial trajectory conditions 
and performs the calculations for a precision-integrated trajectory. The 
input state, options and the tape- or disc-stored quantities are the same 
as in CONW. 
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The targeting option of PINT, REFINE, serves the same basic purpose as 
SEARCH - to try to satisfy end constraints by systematic variation of 
initial conditions. There are some major differences, however. On the 
positive side, REFlNE's sensitivities of constraints to controls are gen- 
erated by integration of variational equations along with the trajectory. 
This method is faster than the secant method used in SEARCH. On the 
negative side, the available control set is limited to the launch param- 
eters listed in Section 1.1.2. The constraint set has been expanded over 
that of the Mark II program to include B*T and B*R or only radius of 
closest approach, but is still very restricted relative to SEARCH. 


Another PINT option enables the user to calculate the sensitivity of the 
trajectory to variations (uncertainties) in constants in the equations 
of motion (e.g. , planetary masses, gravitational harmonic coefficients). 
These sensitivities are calculated by integration of variational equations 
along the nominal trajectory. The same sensitivities are computed, but 
not available for output, in error propagation when treating equation 
of motion parameter uncertainties. Calling out yet another PINT option, 
the state transition matrix may also be calculated by integrating varia- 
tional-equations. The implemented variational equations are described 
fully in Appendix B of Reference 1. 


Some improvements have been made to the set of stopping functions used 
when integrating trajectories in both PINT and SEARCH. The stopping 
functions are computed along with a trajectory and are used to terminate 
the integration. The five functions which are computed (in subroutine 
FSUB) are: 

1. I— i - .03, the Encke rectification criterion 
r 1 


2 . 


r sign (R* V) 


r patch’ • 


the patch -away criterion 
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3* ^^ r patch^’ P atc ^“ to or closest approach criterion 
defined in subroutine PATCH (Ref. 2) 

4. t - t , the time limit criterion 

stop 

5. f ( r stop ) » t * ie criterion for closest approach or given 

radius, from the target body. 

In the Mark II program, the patch-away criterion was not signed, the 
closest -approach function did not exist, and the patch-to and radius- 
stop functions were merely distance differences. The newer formulations 
eliminate ambiguities and permit lengthening the integration step size 
for interplanetary transfer trajectories. 

1.2 ERROR PROPAGATION 

The fundamental questions answered by the Mark IV Error Propagation 
Program are, "Given a space trajectory and a set of measurements of a 
specified type and quality, 

1. How well can the trajectory be determined? (navigation) 

2. How well can measurement biases and equation of motion 
parameter uncertainties be determined? 

3. How do measurement biases and equation of motion parameter 
uncertainties affect the quality of trajectory determination? 

4. What effects do navigation errors have on midcourse guidance 
requirements and how do guidance execution errors affect 
navigation?" 

In order to answer these and other questions, the Mark IV Error Propa- 
gation Program utilizes statistical principles to determine ensemble 
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characteristics of trajectories in the vicinity of the nominal 'trajectory. 
A covariance matrix of state estimation errors is assumed to represent 
these ensemble characteristics. This covariance matrix is manipulated 
by the program according to Schmidt-Kalman filter theory using optimal 
weighting of all measurements. A brief summary of these manipulations 
as implemented in the Mark IV program will be presented in this section. 
More detail and derivations will be found in References 1, 2 and 3. 

The state vector X, consists of 6 cartesian components of position and 
velocity, k equation of motion or dynamic parameters and & measurement 
biases. (The term "measurement biases" includes location errors and 

A, 

time biases.) If X represents the estimate of state and x = X-X repre- 
sents" the error in estimate of state, the covariance matrix of state 
estimation errors, P, is 

P = E (xx T ) 

T 

where E is the statistical expectation operator and X means X-transposed. 
The user must supply P to the program initially. Except for the part of 
P which represents the trajectory estimation error distribution (upper 
left 6x6 matrix) the program expects 'P to be initially diagonal. That 
is, the dynamic and measurement biases are assumed to be uncorrelated 
initially. 

The answers to the questions posed earlier are to be found through inter- 
pretation of the P-matrix. We assume that the mean state error is zero. 
The elements of P represent statistically 




c. 

r 


a . 
J 


where is the standard deviation of state error component x^ about a 

zero mean value and the p . . are coefficients denoting the correlation 

i J 

between x.^ and x ^ . In general, the smaller the standard deviation of a 
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component of state error, the better that state component is estimated. 

The probability that a sample state error vector would lie within the 
hyper-ellipsoid defined by the covariance matrix could be (but is not) 
computed. 

The program computes the trajectory state and its corresponding portion 
of the covariance matrix referred to the mean equator and equinox of 

1950.0 coordinate system. This system is not very meaningful for output, 

so the program supplies simply-computed performance measures for positions, 
RMSP, and for velocity, RMSV, These are defined by 

\ 

RMSP - Kl + P 22 + P 33 “ fil + + a z 2 ■ 

BMSV - K* + P 55 + P 66 ‘fj + CT vy 2 + • 

The i-th diagonal element of P represents the variance of the i-th state 

error component about a zero mean value. At the end of a processing 
interval and in special output the covariance matrix is transformed into 
other coordinate systems than EE5Q and printed out. 

Between measurements or other events, P is propagated according to the 
equation 


P(t ) 
n 


e(t 


n 


; t n-l )P(t n-l> 


- T 

* <t 


t 

n- 


1 


) 


where P(t) means "P at time t" and ®(t : t ,) is the state transition 

n n-1 • 

matrix from time t . to time t . The state transition matrix is par- 

n-1 n 

titioned as shown. 


4(t 


n 


5 ‘n-P 


0 I 

0 0 


0 

0 

I 
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In this partitioning, tp represents the trajectory state transition 

matrix, dimensioned 6X6. It is computed in closed-form as an average 

conic transition matrix between the states X(t n ) and X(t ). For 

n-1 n 

details about the computation of cp , see subroutine PHIZ in Reference 2. 

X 

The 6xk matrix represents the sensitivity of the trajectory state to 
variations in the dynamic bias parameters. It is computed by integrating 
variational equations along the nominal trajectory from t ^ to t^. The 
other elements of 9 are either null matrices or identities. The null 
matrix in the first row represents the fact that the trajectory error is 
insensitive to measurement biases in the interval between measurements. 
The remaining elements of $ represent the assumption that the non- 
trajectory state elements are time -invariant biases. 

At a measurement, P is changed according to the equation 

+ - -T -T -1 - 

P = P - P H (HP H +Q) HP • 

where (+) means ''after processing the measurement" and (-) means "before 
processing the measurement". H is the gradient of the measurement with 
respect to the state at the time of the measurement (see Appendix A of 
Reference 1 for complete derivations of H) and Q is the variance of the 
measurement's random error. It is assumed that the measurement errors 
are uncorrelated. If there are several measurements at the same time, 
these are processed individually so that the inverse indicated in the 
above equation is computed as a scalar reciprocal and the gradient, H, 
is a single -row matrix or vector. 

The Mark IV program contains an option for considering the effects of 
dynamic or measurement biases on the state estimation without including 
them as additional state components. The "consider" option is implemen- 
ted by retaining and manipulating the correlations between the "solved- 
for" state and the "considered" state. If the correlation matrix for 
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the "considered" dynamic biases is denoted C and for measurement biases, 

ux ’ 

C , the manipulation equations between events are: 


P(t ) ** $ P(t ,) & T + * c (t .)U T +UC (t ,) * T +UDU T 

n ' n-1' ux v n-1 ux v n-1' 


C 

ux 


< C n> 


$ c 


ux 


(t ) + UD 
n-i 


C (t ) = 

vx v n' 


i C (t ,) 
vx' n-1' 


In the above equations, U is the sensitivity of the state at t to the 

vector of "considered" dynamic biases over the interval from t , to t . 

J n-1 n 

§ is the state transition matrix from t , to t . D is the (diagonal 

n-1 n 

and constant) covariance matrix of "considered" dynamic biases. .A 
derivation of the above equations may be found in Reference 3. Both C 
and C are considered to be zero initially - indicating that the initial 
state and bias errors are statistically uncorrelated. At a measurement, 
the "consider" option is implemented by the following equations: 

Y = HP”H T + HC" G T + GC" T H T + GWG T + Q 
ux vx 

P* - p“ - (p“h T + c“ G T ) (Y) (p"h T + G ■ g t ) 

vx vx 

C + = C - (p"h T + c" G T )(Y) -1 (HC + GW) 

vx vx vx vx 

C + = C" - (P~H T + c" g t )(y) -1 (hc ) 

VX VX VX ' vx' 

In these equations, G is the gradient of the measurement with respect to 
the "considered" measurement bias vector and W is the (diagonal and 
constant) covariance matrix of considered measurement bias errors. A 


1-13 


PHILCO 


SPACE & RE-ENTRY SYSTEMS DIVISION 
Philco-Ford Corporation 



TR-DA2154 


complete derivation of these equations may be found in Reference 3. The 
above equations for updating P, and in time and changing them 
at measurements are the "heart" of the Mark IV Error Propagation Program. 
Most of the other computations in the program are, relatively speaking, 
"programming details". 


The only event at which P is changed other than the measurement 
event is a simulated midcourse guidance maneuver. It is necessary, in 
considering midcourse guidance, to carry along another 6x6 covariance 
matrix, PAR, in addition to P. PAR represents the distribution of tra- 
jectory errors about the nominal trajectory. These trajectory errors are 
assumed to result from random errors in guiding to the desired nominal 
trajectory. The PAR matrix is updated in time by means of the state 
transition matrix, but naturally does not change at measurements. Both 
P and PAR are changed at a midcourse guidance maneuver simulation. The 
assumption of an impulsive velocity maneuver is made according to one of 
three available guidance laws. This assumption leads to a change only 
in the velocity portion of the P-matrix. Specifically, the P-matrix is 
changed at a midcourse maneuver only by the addition of a 3x3 matrix 
representing the velocity uncertainty due to errors in executing the 
maneuver. 


+ 


= P + 


0 


3x3 


0 


3x3 


°3x3 k2 ‘ E < e * T > 


The factor k defines the accuracy presumed in monitoring the correction. 
The PAR matrix after the maneuver is computed to be the sum of the tra- 
jectory navigation uncertainty prior to the maneuver and the execution 
errors. f~ ~l 


PAR + = 


3x3 


3x3 


°3x3 E < ee > 
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X 

Formulation of E( ee ) is accomplished by transforming resolution, point- 
ing and proportional maneuver errors into the program' s cartesian system. 
The theory and equations for guidance calculations are described fully in 
Reference 1. 

The measurement simulation capabilities of the Mark IV. program include 
observations from Earth-based tracking stations, beacons located on the 
Moon or on planets, and devices on board the spacecraft. Storage dimen- 
sions limit the number of stations to 12 and beacons to 10 for any one 
case. The beacons must all be located on the same body. A list of the 
available measurements will be found in the capability summary which 
concludes this section. Formulae for the various measurements in terms 
of the state and station or beacon locations will be found in Appendix A 
of Reference 1. The station and beacon view times are computed by the 
program from the stored nominal trajectory before proceeding with error 
propagation calculations. During error propagation the stations or 
beacons are assumed to observe (take measurements on) the spacecraft 
when it comes into view and then at regular intervals specified by input 
until the spacecraft is occulted. The simulated onboard measurements 
also occur at regular intervals determined by input. 

Midcourse guidance maneuvers (as many as five) are simulated at pre-set 
times along the trajectory. The available guidance laws are: 

1. fixed time of arrival, 

2. fixed energy at arrival, and 

3. minimum fuel, variable time and energy. 

Another option provides the capability of (linearly) propagating the 
P-matrix (and, if guiding, the PAR -matrix) to the end of the trajectory 
to predict the covariance matrix of miss vector (R*T, B- R, time of 
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arrival, v , X and i) errors. The resulting miss vector errors are 

p 

those which would be realized if no more observations (or, for guidance, 
no more correction maneuvers) were made. The expected midcourse velocity 
correction requirement and direction of the critical plane normal vector 
are provided at each output point when the guidance option is selected. 

The theory for these computations is found in Reference 1. 

The program's storage requirements set the P-matrix dimension at 900 cells. 
This permits the treatment of a 30-element state vector for any one case 
(30 x 30 =900). The number of dynamic or measurement biases considered 
but not solved-for may be 100 or less. The principal limitation on the 
number of state components imposed by the 900-element P-matrix is 

(number of solved-for) x (number of solved for + number considered) £ 900 

The special output capability (SPOUT) performs calculations not normally 
required or performed in error propagation. A special output file is 
written on tape or disc by option in ERP. This file contains the tra- 
jectory history and the P-matrix history relative to cartesian EE50 
coordinates, as well as RMSP and RMSV. When SPOUT is called, the file 
is read in and transformed according to input option specification for 
selective output. RMSP and RMSV may be printer-plotted versus time as 
may be any of the solved-for bias standard deviations. A case and record 
number connects the SPOUT file with the ERP output and permits the user 
to be selective in his SPOUT calculations. 

The following summary provides a short but comprehensive list of the 
capabilities of the Mark IV Error Propagation Program. 
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CAPABILITY SUMMARY 


I Trajectory Generation 
A. Start-Up 

1. Determines interplanetary transfer trajectory 

a. massless planet solution 

b. matched-conic finite influence solution 

2. Determines necessary injection conditions 

a. circular parking orbit assumed 

b. computes launch control parameters (see B.l. below) 

c. computes cartesian and spherical injection state 


B. Search 

1. Controls 

a. Cartesian state (x, y, z, vx, vy, vz) 

b. spherical state (r, lat, Ion, v, y, az) 

c. launch control parameters (insertion azimuth, insertion 
time, park time, Av azimuth, Av elevation, Av magnitude) 

2. Constraints 

a. end cartesian state (6 components) 

b. end spherical state (6 components) 

c. miss vector (B*T, B*R) 

d. orbital elements (v , r , i, Q, U) ) 

co p p 

e. time of flight 

3. Search run options 

a. steepest descent constraint satisfaction 

b. secant method gradient with hold-fixed option 

c. patched-conic trajectory model 

d. perturbed, integrated trajectory model (see D.4. below) 

e. midcourse targeting 

f. multiple control state coordinate -systems 
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X.B.3 (continued) 

g. multiple constraint coordinate systems 

h. automatic constraint scan in one control 

C, Patched-Conic Trajectory-Write (CONW) 

1. Generates and writes interpolation coefficients on tape or 
disc file 

2. Writes critical records of patch-points and end points 

3. Generates and writes end point miss sensitivity matrix 

D. Precision Trajectory Integration (PINT) 

1. Integrates perturbed trajectory 

a. writes interpolation coefficients (see C.I., 2., 3. above) 

b. integrates variational equations for state transition matrix 

2. Refines initial conditions 

a. steepest descent constraint satisfaction 

b. integrates variational equations for gradient 

c. launch parameter controls (see B.L above) 

d. constraint options (B*T, B*R, r , v , time of flight, 

p CO 

target vector, earth-return point, minimum injection Av) 

e. target body orbital constraint reference system 

3. Computes motion parameter sensitivities 

a. interpolates trajectory from tape 

b, integrates variational equations for sensitivities 

4. Trajectory calculations 
Encke's Method of trajectory calculation 
Adams’ fourth-order method of numerical integration 
interpolated stopping functions (Encke rectification, 
patch-away, patch-to, closest approach, stopping radius 
time of flight) 
perturbations considered: 

4 zonal harmonics of Earth's gravity 

5 longitudinal harmonics of Earth's gravity 
7-body gravitational attraction 
lunar triaxiality 
Earth's atmospheric drag 
solar radiation pressure 
tangential thrust 
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I.D.4 (continued) 

e. center-shift at fixed spheres of influence 
E. General 

1. Optional lunar and planetary ephemeris 

a. approximate, mean-element package 

b. JPL ephemeris tape interpolation package 

2. Units of computation 

a. time in seconds 

b. positions in kilometers 

c. velocities in kilometers /second 

3. Coordinate system for computation 

a. Earth's mean equator and equinox of 1950.0 

b. cartesian 

II Error Propagation 

A. Measurements (Random Error Sources) 

Earth-based tracking stations (up to 12) 

a, range 

b, azimuth and elevation 

c, right ascension and declination 

d. Minitrack (direction cosines) 

e . range rate 

f, azimuth and elevation rates 

g. right ascension and declination rates 

h. direction cosine rates 
Moon- or planet-based beacons (up to 10) 

a. range 

b. range rate 

c. azimuth and elevation from the vehicle 
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II .A (continued) 

3. On-board star and planet measurements 

a. height 

b. height rate 

c. planet subtended angle 

d. latitude and longitude (from inertial platform) 

e. sextant (star-planet angles) 

B. Deterministic Error Sources 

1. Equation of Motion Parameters 

a. astronomical unit conversion 

b. lunar and planetary masses 

c. zonal harmonics of Earth's gravity 

d. longitudinal harmonics of Earth's gravity 

e. harmonics of the Moon's gravity 

f. Earth's atmospheric drag coefficients 

g. solar radiation pressure coefficient 

h. venting thrust magnitude 

2. Deterministic Error Sources in Measurement 

a. tracking station measurement biases 

b. tracking station location errors 

c. tracking station time bias 

d. beacon measurement biases 

e. beacon location errors 

f. vehicle time bias 

g. onboard height and height rate biases 

h. velocity of light uncertainty 

C. Program Implementation 

1. Minimum variance estimation technique 

2. Linear propagation of the covariance matrix 

a. closed-form conic transition matrix 

b. computed in inertial Cartesian coordinates 
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II. C (continued) 

3. Automatic search for station on-off times 

a. occuLfcations 

b. artificial horizons and zenith limits 

4. Nominal trajectory interpolated from tape 

a. patched conic 

b. precision integrated 

5. Extra output tape 

a. normal output of rms position and velocity without tape 

b. covariance matrix output 

c. state output various coordinate systems 

d. automatic plotting of specific parameters 

D. Treatment of Deterministic Errors 

1. Treat as though solving for the error 
a. 900-element covariance matrix 

2. Treat as though only considering the error's influence 
a. consider up to 100 error sources 

3. Equation of motion error sources 

a. sensitivities obtained by integrating variational 
equations about nominal trajectory 

E. Prediction and Guidance 

1. Prediction to the end point 

a. uncertainty in miss parameters due to navigation errors (P) 

b. uncertainty in miss parameters due to trajectory errors (PAR) 

2. Midcourse guidance 
a. Guidance laws: 

fixed time of arrival 
constant target-relative energy 
minimum velocity correction 
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II. E. 


PHILCO 


(continued) 

b. Correction errors: 
pointing 
resolution 
proportional 
monitoring 
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SECTION 2 
PROGRAM USAGE 

The Mark IV Error Propagation Program consists of six related sub -programs 
which may be run individually or serially. These sub-programs are: 

1. PLANET (START-UP) for generating approximate interplanetary 

trajectory information and initial conditions 

2. SEARCH for improving initial conditions in order to satisfy 

specified end constraints 

3. CONW for writing an ephemeris file of the vehicle's patched- 

conic trajectory 

4. PINT for writing an ephemeris file of the vehicle's integrated 

trajectory 

\ 

5. ERP for performing error propagation calculations 

6. SPOUT for calculating and displaying selected quantities 

after error propagation. 

The organization of the program by subroutine is shown in Figure 2-1. 

The purpose of each subroutine is listed in Appendix E. Descriptions 
of the subroutines are to be found in Reference 2 as amended under this 
contract. 

Each sub -program has need of planetary ephemeris information. This 
information may be supplied to the program by means of a self-contained, 
mean-element ephemeris (subroutines ANTRI, EXPAND, EL2EX and UPDATE as 
shown in Figure 2-1) of by means of JPL's ephemeris tape (interpolated 
with subroutines ANTRI, DEPHEM , BUFF IL and INTCOF). The taped ephemeris 
slows program execution and requires more storage, but is more precise 
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chan the self-contained ephemeris. Replacement of decks is required 
in order to change ephemeris options. 


2.1 INPUT DATA 

A data run consists of one or more cases submitted to computer operations 
with the Mark IV program deck for the purpose of obtaining printed output 
information computed by the program. Each case consists of a call of one 
of the six main sub -programs. The first card of any case contains a 
number (punched into column 5) which specifies which of the main sub- 
programs is to be called for that case. The code for this card is: 

0 end of run 

1 PINT 

2 ERP 

3 SPOUT 

4 CONW 

5 PLANET 

6 SEARCH 

Several precautions must be taken in ordering cases. First, a PINT or 
CONt 1 / case must be executed sometime prior to an ERP case so that the 
vehicle's ephemeris file required by ERP has been written. Second, a 
tape or disc file must be written by ERP before SPOUT can be run. These 
files may be written in a previous run if stored on tape. In this event, 
the user must be sure to have the right tape mounted for the later run. 
Another precaution is to avoid executing PINT, PLANET or SEARCH cases 
after ERP, SPOUT or CONW cases in any one run. The first three sub- 
programs use WCOM as an input array and the latter three use WCOM as 
a working array. Incorrect usage could erase expected input data, 
causing strange results. 

Case data follow the first card for each case. These data replace or 
augment standard case data stored in the program by BLOCK DATA routines. 
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The stored data values may be found in computer listings of INPCOM and 
WCOM or in Reference 2. Data contained in the program at the beginning 
of a run will remain and be used until changed by input. In other words, 
if a great many input quantities are specified in case no. 1 and desired 
for other cases in the run, they need not be specified again for later 
cases in the run. When changed for any case, the changed data values 
will remain in use for the remainder of the run or until changed again. 
The only exception to this philosophy is in ERP. 

Input data are read into the two COMMON arrays, INPCOM and WCOM. Gener- 
ally, the user must supply the address or location of the array to be 
filled and the numerical value to fill that address or location. Input 
(Ascription for the Mark IV program consists primarily of lists relating 
locations to parameters required by the program. Appendix D of this 
manual is a summary of input requirements of the modified Mark II pro- 
gram’s input requirements and sample input/output may be found in 
Section 5 of Reference 1. The remainder of the present section will 
be given to input and output examples of the modifications: the start-up 

and search capabilities. 

2.1.1 INPUT FOR START-UP 

The start-up option is called with "5" in column 5 of the first card. 

Case data are then read into the WCOM array through subroutine ROVLEY. 

An input example is shown in Figure 2-2. The launch and target planets 
are specified by fixed-point inputs and the body code of Table 1. 

IN(1) launch body number (1 for Earth) 

IN (2) target body number (5 for' Mars) 

The remaining case data are input by locations and corresponding numer- 
ical values in a 4(13, E12. 8) format. (See Table 2.) These are 
primarily: 
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departure date (7307.24, 0. means 0* 1 , July 24, 1973) 
arrival date (7402.16, 0. means 0 , February 16, 1974) 

although it is also necessary to have specified the number of matched- 
conic iterations to be performed (X(405) = 1.0) and the desired miss 
vector components at the target (B*T = X(373) =-4930.0 km and B*R = 

X(374) = 7210. km). Parking orbit parameters at launch must also be 
specified. 

insertion latitude (X(121) = 28.5°) 
insertion longitude (X(122) =-80.5°) 
insertion azimuth (X(123) = 67.74°) 
insertion altitude (X(124) = 184. km) 

These parameters enable the program to compute specific injection con- 
ditions which are meaningful relative to the parking orbit assumption. 

An output sample resulting from the input example of Figure 2-2 is shown 
in Figure 2-3 and described in section 2.2.1, 

2.1.2 INPUT FOR SEARCH - 

The search option is called with "6” in column 5 of the first card. Case 
data are then read into the WC0M array through subroutine R0VLEY. The 
principal data required for a search case are: 

initial date and time (locations 99 and 100 or 101 and 102) 

-initial trajectory state (311-316, 317-322, or 323-328) 
control selector (310) 

control limit levels (331-336, 337-342, or 343-348) 
constraint selectors (351-356) 
desired constraint values (361-380) 
constraint error tolerances (381-400) 
option selectors (scattered locations) 
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The complete list of input data required by SEARCH is found in Table 3. 


An input example is shown in Figure 2-2. The user must supply the launch 
and target body numbers unless these have been supplied to the program in 
an earlier TINT, PIANET or SEARCH case of the same run. The input example 
for PLANET shows how to set IN(1) = 1 and IN(2) — 5 which tells the program 
that the launch body is Earth and the target is Mars. The SEARCH case 
shown was run after a PLANET case, so the launch and target bodies were not 
re-specified. The next line (card) shows locations 111, 112 and 113. These 
represent (respectively) the time (0. DH.MS) from initial date at which the 
variation or search is to begin, the precautionary stop time (21000. DH.MS 
= 210 days), and the precautionary distance (3400. km) from the target at 
which the trajectory must stop if neither stop time nor closest approach has 
been reached. The next card contains the fraction of control limits to use 
in generating partials (X(146) = .0001) and an upper constraint tolerance 
factor (X(175) ** 100.) used in iterative testing. The next two cards 
specify eight options. 

(301) 10. means that the iteration is to be terminated after 10 
iterations even if convergence has not been achieved 

(302) 2. means that the program is to compute its own factors 
for scaling the gradient 

(303) 1. means that the gradient is to be computed by finite 
differences in FNDMXN rather than being supplied 

(304) 2. means that the iterative steps are to be controlled by 
numerical estimation of curvature when minimizing 

(305) 0. asks for no extra output 

(306) 2. means that the. gradient is to be used for two more 
iterations in addition to the one in which it was computed 

(307) -1. cancels extra output of initial conditions at each step 

(308) 0. asks for the patched-conic trajectory model 
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The next card sets 3. into location 309 and 11001. into location 310. This 
3. means that the control set is to be the launch parameter set for which 
initial conditions are located in 323-238. These initial conditions are 
omitted from input because they have been set by the earlier PLANET case. 

Other parameters which would have to be input for SEARCH if they were not 
set in the PLANET case are the park orbit insertion date (locations 99 and 
100) and the park orbit specifications (locations 121-124). The 11001. in 
location 310 specifies the control variables to be launch time, park time 
and injection velocity impulse. The next card specifies the control limit 
levels on launch time (100. seconds), park time (5. seconds) and injection 
velocity impulse (.01 km/sec). The next four cards specify the constraint 
set, desired values and tolerances. 

Selection Desired Value Tolerance 

(351) 13. for B* T (373) B*T is -4930. km (393) B'T is 10. km 

(352) 14. for B-R (374) B*R is 7200. km (394) B*R is 10. km 

(353) 20. for flight time (380)t f is 206 d 14 h 30™ 0 s (400) t f is 100. seconds 

(354) 0. ends the constraint set length at three 

The next card selects the search option (S(418) = 0.) as opposed to the 
scanning option and selects the target-fixed constraint coordinate reference 
system (X(419) =3.). A blank card ends the input and leads to execution 
of the case. An output example resulting from this input example is shown 
in Figure 2-3 and described in Section 2.2.2. 

2 . 2 PROGRAM OUTPUT 

Print-out from the Mark IV Error Propagation Program is basically the same 
as that of the Mark II prog-ram except, of course, for the addtion of start- 
up and search output. Examples of the Mark II program's print-out are given 
and explained in Section 5 of Reference 1. Examples of the start-up and 
search print-outs will now be given. 

2.2.1 OUTPUT FOR THE START-UP CAPABILITY 

The start-up capability's output is computed and printed primarily in 
PLANET. Figures 2-3 and 2-4 show this print-out. The first block of print- 
out for any case shows the data input for that case (see 2.1.1) and is labeled 
"OVERLAY INPUT." The massless planet solution is then headed "APPROXIMATE 
TRAJECTORY DATA." Later iterations are headed "MATCHED- CONIC ITERATIVE 
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3. HAT 

k rt/a 

HYP EXCESS SPEED 

b.874 

KM/S 

CJ IENERGY) 

14. 64 

K2/3? 

CJ ( energy > 

«.l 6 

K?/S2 

latitude* Ecliptic 

2 >.«56 

nE^»a 

LATITUDE* ECLIPTIC 

-?R.5S6 

OEGS 

LATITUOE. ORuITAr 

23. 154 

C-EOb 

LATITUDE* ORBITAL 

-28 2 1 3 

DEG5 

DECLINATION 

35.414 


0LCLIN4TI0N 

- 2.5*1 

DFBS 

LONGITUDE* EiLtEoU 

16.777 

nE^6 

LONGITUDE* ECLfEQU 

i9. ns 

OFGb 

LONGITUDE* QKH*RaD 

92.97*1 

nEOb 

LONGITUDE* ORB * RAD 

-53 3 9 

DEGS 

RIGHT ASCENSION 

2S,3i 

l)EGb 

RlUMT ASCENSION 

45.7 6 

decs 

SUN-ASYMPTOTt ANR. 

87. 7«- 

r,£Gb 

SUN-aSYmPTOTE ANG, 

1 '1.771 

DECS 



SPECIFIC 

LAUNCH 0 A 1 A 




INSERTION. DATE AND par* ORHIT DtSCRTPTION 

JUL 2*. 973* R rtHS. 36 MIN* 43.039 SEC * JULIAN DATE ?44}887 .65884) 89 

ALT 1 *84000000t *02 L.A7 2,8500(JOQ0r.(il ^ON-H . f}5'>QOQOOF*0 1 VA^ 6. /7400000E*01 (N<- J .55 70l4£»5E*fil - <*HA 7. llo7732Zfc-*Ol 


INJECTION CONDITIONS At EaPTH 
JUL 24* :973t V HRs, 36 MIN, 
DLA^ 0 . "TL 


X-2.2t'352443E*o3 Y-S.2447b22iP-.i3 2-3.?7i'94234E*03 

R 6.56216500E*01 OEC-2.yR97H199r.oi RA-1 . 12 '«9l 7lE*n2 

SMA-2, 83*2901 7E*04 £LC 1 . 23 152764F*'>n INC 3 . 5465ft 7 13E*0 1 
C3 1 .4C-6359H9E*0l THET-9 , 432D4724F-09 PERV 1 , 1 6425457E ♦ 0 1 


S6. -‘23 SEC JULIAN DATE ?*4 1887. 90064957 

PfKT 3.6l^lR369E*03 DVAZ-5 .8153599 IE-U DVEL-^.l 1229914E-13 r»£LV 3«84378950E.nQ 


DA 9.31418175E+00 DY-6,0696n 156 f *00 02 3. 4575 7520E.00 

V 1 ,164254S7E*01 PTH-5.2 o 50 ifllBF-09 A2 6, 9966591 4E«4l 

L*N 3.0 1 024531E+02 APF 3.00786337E*02 RCA 6,5621 6500L*03 
SLR 1,46436525E*04 JMPV 3.84878950F* DO TPER- 1 . 07391222E-1 2 



Figure 2-4. OUTPUT FOR START-UP 
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SOLUTION." The next two lines show the departure and arrival dates and 
planetary positions at those dates. The coordinate system for position 
output here is the mean ecliptic and equinox of launch date. Radii are 
printed in a.u, and latitudes and longitudes in degrees. Heliocentric 
orbit data are printed next to describe the heliocentric transfer trajectory. 
Transfer time and transfer angle refer to the interval between departure 
patch and arrival patch - and the sphere of influence has zero radius for 
the massless planet solution. The only other explanations needed for this 
block are that "TRUE ANOMALY" refers to the heliocentric orbit at departure 
patch and that the ascending node, inclination and argument of perihelion 
are referred to the mean ecliptic and equinox of launch date. The "ASYMP- 
TOTIC DATA" describe the hyperbolic excess velocities at departure (LAUNCH) 
and arrival. Declination and right ascension of the asymptotic (hyperbolic 
excess velocity vector) are referred to the Earth’s mean equator and equinox 
of launch date. The only other explanation needed is that the sun-asymptotic 
angle is the angle between the p lanet-to-sun line and the hyperbolic excess 
velocity vector. 

"SPECIFIC LAUNCH DATA" describe the park orbit, insertion and injection 
conditions which connect the parking orbit with the departure hyperbola. 
Insertion date is the time of park orbit initiation. Injection conditions 
refer to the initial conditions for the departure hyperbola. The reference 
frame is EE50 for all cooordinatized data. The following list explains the 
remaining symbols. 

Parking orbit altitude (km) 

Parking orbit insertion latitude (deg) 

Parking orbit insertion longitude (deg) 

Parking orbit insertion velocity azimuth (deg) 

Parking orbit inclination (equator of date, deg) 

Greenwich hour angle at insertion date 

Incremental insertion azimuth control (rad) 

Incremental launch time (sec) measured from insertion date 
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ALT 

LAT 

LON 

VAZ 

INC 

GHA 

DLAZ 

DTL 
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PRKT 

DVA2 

DVEL 

DELV 


Time in parking orbit (sec) 

Azimuth of the injection impulse measured CCW from the 
parking orbit at injection (rad) 

Elevation of the injection impulse measured up from 
local horizontal at -injection (rad) 

Magnitude of the injection velocity impulse (km/sec) 


X,Y,Z Cartesian position components at injection (km) 

DX,DY,DZ Cartesian velocity components at injection (km/sec) 


R,DEC,RA Spherical position components at injection (km, deg) 
V,PTH,AZ Spherical velocity components at injection (km/sec, deg) 


SMA 

ECC 

INC 

.LAN 

APF 

RCA 


Semi-major axis (km) 

Eccentricity (no units) 

Inclination (deg) 

Longitude of the ascending node (deg) 
Argument of perifocus (deg) 

Radius of closest approach (km) 


C3 

THET 

PERV 

SLR 

IMPV 

TPER 


2 2 

Vis -viva energy (km /sec ) 

True anomaly at injection (deg) 

Periapsis velocity (km/sec) 

Semi- latus rectum (km) 

Impulsive velocity required to circularize (km/sec) 
Time after periapsis passage (days) 


2.2.2 OUTPUT FOR THE SEARCH CAPABILITY 


The standard output for the search capability is shown in Figure 2-5. 
Input case data are printed in the block labeled "OVERLAY INPUT." The 
next few lines are printed in SETUP and describe the control-constraint 
parameters in effect. For the case shown, the control variables are time 
of launch (DTL) , park time (PRKT) and velocity impulse (DELV). Initial 
values and limit steps corresponding to the control variables are shown 
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in 

1 .C00O0000E-04 

lid 

overlay input 
2.IOOOOOO0E*p4 

113 

3.40000000E*03 

-0 

— u . 

I'»6 

l7o 

l .nooooonoF^o? 

-0 

“0. 

-0 

-o. 

3 'll 

l.(]00uo00flt*01 

302 

7.ooooonoPE*oo 

303 

1.000u0000£*00 

J04 

2.00^0n000c*00 

305 

• 

3(\b 

? , OOflOOO O0F*OU 

307 

-1 .OOOOOOOOf *00 

30^ 

c ♦ 

J09 

3. 000000001- *00 

31P 

1.10010000^*04 

-0 

-o. 

-P 

-u. 

J44 

1 • 00 00 000 0E* 02 

34b 

S.0DC0000PF*0U 

348 

1 .000U0000E-02 

-0 

-0. 

3Si 

1 . 30OOOO 00 1£*0 1 

373 

-4.93o000one-f J 

3V3 

1.000O0nooE*oi 

-n 

-V, 

352 

‘ ,40000000E*A1 

J7h 

7 *?0000000£ 4, i3 

394 

I .OOOOOOOOE*Ol 

-0 

-v. 

353 

2.00000000E*01 


?« 06143000 t*r4 

400 

l.OOOO0i)Ol)E*0? 

-0 

“0. 

354 

• 

~ 

-A . 

-0 

-0. 

-0 

—0 . 

-18 

* 

.w 

3. TOOOOOOOr *00 

-0 

“0 . 

-o 

-0. 


TaWGET THk T* JhCTORY FROM fcAHTrt TO MARS 


CONTROL PARiiMtTFBS 

INITIAL VlUJtb 

L’TL 

, 

PKKT 

'.61218369k *0 * 

l)tL 

i.R487895.»E*0i’ 


limit LEVELS 
l .uoooooooe*o 2 
^.vooonooor+oo 

l.COOOf'OOOt-n2 


CONS TKM MTS 


OtSIHEO V Hits 


TOLER NCFS 


H|» I 

i>rw 


:fl 


~4.930f1000AF*()J 
/. 0000000 k *03 
2 o 6 i 4300 ')t*Q 4 


1 .O 00000 o 0 fc* 0 l 
1 ,C UQ 0000 QL* D 1 
1.6OOOOOf)Ot*02 


JUL Z4» 973. 

JUu 24. *973. 

DLA2 r. 

X-2.2‘'3b2S‘i0t*C3 
K 6.562l650ut-<iJ 
SMA-? • H34P9Q1 7t* '4 
CJ UAJ'^JSyHVfc^ol 


H M«s. 36 MIN, 43 .o3« SEC 
9 Hris» 36 MIN, 5b.. 23 S' C 

DTL PHUT 3.M2U169E»'3 

Y-S. 244701 72r*p3 Z-3.?7u ,42l3F*nj 

PfcC-2.9BV7oi 78 f* 1 oa-I 
F.CC 1 • ?315« > 764 f*“ 0 INC 3 ,S406i,b66F *<»l 
HtT-9.<.33?337ir»»9 P£py 1 . 1 64254S7E * rt l 


Julian date ?*4] 887.usyu*i«9 

JllLlAM DAT- ?44l8b/. 90064957 
DvAZ-b,8lb3599lE-l4 0VFL-<:.-l2299l4F-13 OELV J.8487b950fc*o0 
Dt 9.3141 ,07SE*00 [>Y-b,t>6-6<>3bOE*O0 0/ 3.45757431 t*O0 

V I«16425457E*P1 PTM-b.20S49293F-n9 A/ ©.996609 72fc .01 

l AN J.ql0?4S73F*02 4PF J.nft78A33?F*ft? RCa b .56^ l 6b(>0t * 03 
SL6 1 .46* J652SE*04 ImPV J.8487H9b0F*O0 TPER-1 .0 740*73l£-12 


ITEM. no. 

0 

OnTpuj S 
CONSTRAINTS 

control inc 

0.00000000 

-517.2982 

H6.220r,2*nj 

36l2,18J'>9'*54 
-?154V.a4?l 
-4.P27' sfl39 

3.848*895 
6643 • J i 2b 
-.00032891 

IIEH, NO. 

1 

CQnTroi b 
Constraints 
CONTROL INC 

86.22082403 

39.9275 

-21.67043382 

3607.956 'iol5 
“29 .7345 
l. .71 . 159 

3.84846>&9 

847.5890 

.00002387 

i TEH, hO. 

2 

controls 
constraints 
control INC 

64.549R9 21 
-5,3-40 
.101rt50l3 

3609. 2*05174 
13!>*37 9 
-.006 0516 

3.848484*6 
-61 ,bb28 
.00000115 

ITEK. «0. 

3 

CDnTrulS 
Co *STRAINTS 

64.6bA94f 34 
-.314"* 

36-9 • ' 2225668 
-.7563 

3.84848601 

-9,089b 



«* IT CO VERGES IN 1 IIERAIIONS ***** 
INITIAL CONDITIONS FARTm -CENTERED 


JUL 24* 1 9 73* 9 m«S. 37 MIN, -7.612 SEC 

OLAZ ... DTL 6. 465094 03F*d! PRUT 3, 60*i ?226E*03 

*-?•! 9fl0323RE*f,3 Y-S.2422©713 F*a3 2-3 .? 7«2930I>F ♦ 03 

R 6 *5821650 Ut*o3 I)tt-2.9971d804r*01 SA-1 . 127b2524E*f 2 
SMA-2.8ib7l6V6t*^4 EC C 1 .23141 1I4F*O0 INC 3 .546© i 88 oE* ol 
C3 1 .41 56522«t*0l THf-:T-9.55b8b7*lR-i)9 PERV ] • 1 64Z?4 J RE*GI 


Julian date 2441887,90136125 
DvA,-b.8Ib3599lE-14 DVEL-^. *>1 229914F-1 3 OELV 3. 84848561 L*0O 
0* 9 ,3p7t)6807E*OO DY-<>.06P469S3E*00 02 3.435604925*00 

: V 1 ,16*22M8E*0I PTH-S,27326028E-f)9 A7 7 . 0083 J003E *0 1 

IAN 3.0 1 Z95380E* 02 APF J.0057rt289F*l>2 RCA 6.56216500^*03 
SLR 1 •‘♦64 28881 E* 04 IMPV -J. 84R4R561F* DO TPER-1 ,0880J73fit-lZ 


TEQMIajAL CONDITIONS MaHS -CENTERED 


Ftrt Ihi 1974. * 

x 4*776. i9S7t*03 
K 4^9' H4884oE*o3 
SMA-5*303C0326t*Cl 
C3 H.l' 44?647t*00 


H«S, 8 MIN. 6.*r 2 SEC 

Y-7,75«46309f*02 l 8 . 254m 27E*02 
DtC 9.67838827 f*<io »A-9 .22o©7224£*no 
ECC 1 .92b6. , S40F*(10 INC 2. 33bl696lF*0l 
ThLI- 7 , 64854j87r-i J PERV 5.{)6l23004F*00 


JULIAN DATE 7442094. 50563312 
l>- 1.06165184E*00 DY 4, 6030H009F + 00 D/-1 . 8 167S496E* 00 

V 5.06l23004E*00 PTH-4,5474735lr-l3 AZ 1 . 11 354544b* 02 

IAN 1.94040315E*02 aPF 1 .5490i842E*OH RCA 4, 9084«B48E*03 
SLR I.43603004t*04 ImPV ^.3022n798 F *00 TPEP-l , 49842142^-16 


Figure 2-5. OUTPUT M)R SEARCH 
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in their appropriate units, i.e., seconds, seconds and kilometers/second 

A A 

respectively. The constraints for this case are B * T, B • R and time of 

A A 

flight. The desired values and tolerances for B *T and B ' R are in units 
of kilometers, while the desired value of flight time is printed out in the 
(days) (hours) . (minutes) (seconds) format used for input.. The flight time 
constraint tolerance is in units of seconds. Only flight time of the 20 
constraints has mixed units for desired value and tolerance. 

4 

The insertion date is next printed if launch control parameters are selected, 
followed by the initial estimate of injection date ana v the initial values of 
the entire launch control parameter set (see 2.2.1 for explanation of sym- 
bols). If either cartesian or spherical controls had been selected, the 
insertion date would not have appeared and the starting date would be 
printed in the calendar format used for input. The input initial state 
would be printed out next. In either case, initial cartesian, spherical and 
orbital components of the trajectory are next printed out in EE50 coordinates. 
(See 2.2.1 for symbol definition and units.) 

A number of blocks are next printed to show the iteration history. The 
iteration number is followed by the control values used for that .iteration. 

The order and meaning of these values follows that of the earlier print-out 
of control variables in effect, i.e. DTL, PRKT and DELV. The constraint 
errors are printed next in the earlier-stated order. The error is defined 
as "desired value minus current value." Units of the constraint errors are 
the same as input units except for flight time, where seconds are seen. The 
third line in the block shows the control increments as calculated by the 
program to improve the constraint error. These increments are in the same 
order and units as the control values. When a control increment is the 
same size as the limit level, a limited control step will be taken. The 
iteration is considered to be convergent when t-he constraint errors are 
each smaller than the corresponding tolerances. 

The final two blocks of print-out show (1) the initial date and initial 
state which correspond to the solution control set, and (2) the final date 
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and final state which were determined by the solution control set and the 
trajectory model in effect. The coordinate system for both of these blocks 
is EE50. The symbols and units for these blocks are identified in 2.2.1 
although the comment "at injection" no longer necessarily applies. 

The user may obtain extra output by setting the key in location 305 non- 
zero. In this case the time and state conditions at trajectory initiation 
and at closest approach to the target body are printed out. The latter 
are printed in the selected constraint coordinate frame at the end of each 
trajectory calculation. Another line then prints out the computed values 
of the quantities listed below. 


BDT 

B * T (km) 


BDR 

B •£ (km) 


TFL 

Time of flight (sec) 


DECS 

Latitude of the arrival asymptote 

(deg) 

RAS 

Longitude of the arrival asymptote 

(deg) 

RAT 

, A 

Longitude of the T-vector (deg) 



The extra output also includes the scaled gradient. 



and the factors by which it is scaled. The elements, of the gradient may be 
obtained by dividing the scale factor into each element of the column to 
which it belongs. 


Extra output may also be obtained for integrated trajectories by setting 
the parameters in locations 115 and/or 116 non-zero. This extra output 
is primarily used for de-bugging purposes. 
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TABLE 1 

CELESTIAL BODY NUMBER CODE 

1 Earth 

2 Moon 

3 Sun 

4 Venus 

5 Mars 

6 Saturn 

7 Jupiter 
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TABLE 2 

INPUTS FOR START-UP 


Case Card: Put 5 in column 5, anything or nothing thereafter 

ROVLEY cards: Integer or real data ended by a blank card 


Location 

Name 


Units 

IN(1) 

IB1 

Launch body number (see Table 1) 

None 

2 

IB2 

Target body number (see Table 1) 

None 

X(121) 

B0LAT 

Park orbit insertion - latitude 

Degrees 

122 

B0L0N 

Park orbit insertion longitude 

Degrees 

123 

B0VAZ 

Park orbit insertion azimuth 

Degrees 

124 

B0ALT 

Park orbit insertion altitude 

Kilometers 

373 

BT 

A 

Desired miss vector component, B'T 

Kilometers 

374 

BR 

/> 

Desired miss vector component, B*R 

Kilometers 

401 

DL1 

Launch date (year, month, day) 

YM.D 

402 

DL2 

Launch date (hour, minute, second) 

HM.S 

403 

DAI 

Arrival date (year, month, day) 

YM.D 

404 

DA 2 

Arrival date (hour, minute, second) 

HM.S 

405 

YET 

Number of matched-conic iterations 

None 

412 

DIAZ 

Incremental insertion azimuth 

Radians 
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Table 3 

INPUTS FOR SEARCH 

Case card: Put 6 in column 5, anything or nothing thereafter 

ROVLEY cards : Integer or real data ended by a blank card 


Location 

Name 

Meaning 

Units 

IN(1) 

IB1 

Launch body number (see Table 1) 

none 

2 

IB 2 

Target body number (see Table 1) 

none 

X(101) 

YW 

Initial date (year, month, day) 

YM.D 

102 

YF 

Initial date (hour, minute, second) 

HM.S 

111 

TSECI 

Time from input date above at which 
the search is to originate (used 
only for cartesian or spherical 
searches) 

DH.MS 

112 

TSTP 

Time at which the trajectory is 
to stop even if the target body has 
not yet been encountered 

DH.MS 

113 

RSTP 

Distance from the target body at 
which the trajectory computation ends 
even if closest approach or stop time 
have not yet occurred 

kilometers 

146 

ESX 

Fraction of control limits to be used 
in generation of partials 

none 

175 

EPSM 

Upper constraint tolerance factor 

none 

301 

TRIES 

Maximum number of iterations before 
giving up the search or the number 
of scanning steps 

none 

302 

SC0PT 

Scaling option key for gradient 
computation (use 2. for automatic 
scaling) 

none 

303 

GROPT 

Gradient option key (use 1. for 
automatic differencing) 

none 

304 

COPT 

Convergence control option key 
(use 2.) 

none 
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Table 3 (continued) 


Location 

Name 

Meaning 

Units 

305 

XTRA 

Extra output key for gradient and 
target periapsis conditions (0. to 
suppress, 1. to obtain output) 

none 

306 

SAVG 

Number of trials to use the gradient 
after the trial on which it is 
computed 

none 

307 

XWOUT 

Extra output key for initial conditions 
(0. for extra, -1. for none) 

none 

308 

HOW 

Trajectory option key (0. for patched- 
conic model, 1. for perturbed, 
integrated) 

none 

309 

TYPE 

Control set option key 
(1. cartesian, 2. spherical, 
3. launch parameters) 

none 

310 

XNI 

6 ' -1 

Control selector key (10) 1 

where is assigned by the following 

chart : 

none 


N. 

i 

TYPE = 1. 

TYPE = 2. 

TYPE = 3. 

N 6 

X 

r 

DLAZ 

N 5 

Y 

lat 

DTL 

A 

Z 

Ion 

PRKT 

N 3 

VX 

V 

DVEL 

N 2 

VY 


DVAZ 

N 1 

VZ 

az 

DELV 


N. = 0. omits the i-th variable from the control set and 

i 

■= 1. includes it. For example, 11001. selects the 
second, third and sixth controls.) 
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Table 3 (continued) 





IF TYPE = 1. 



Location 

Name 

Meaning 


Units 

311 

XIN 

Initial cartesian state values 


kilometers 

312 




kilometers 

313 




kilometers 

314 




km/ sec 

315 




km/ sec 

316 




km/ sec 

331 

XLIM 

Cartesian control limit levels 

(see XNI) 

kilometers 

332 




kilometers 

333 




kilometers 

334 




km/ sec 

335 




km/ sec 

336 




km/ sec 



IF TYPE = 2. 



317 

XIN 

Initial spherical state values 


kilometers 

318 




degrees 

319 




degrees 

320 




km/ sec 

321 




degrees 

322 




degrees 

337 

XLIM 

Spherical control limit levels 

(see XNI) 

kilometers 

338 




degrees 

339 




degrees 
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Table 3 (continued) 


Location 

Name 

Meaning 



Units 

340 





km/ sec 

341 





degrees 

342 



IF TYPE = 

= 3. ~ 

degrees 

99 

YWIN 

Park orbit 

insertion 

date (year, 

month, day) YM.D 

100 

YFIN 

Park orbit 
(YW, YF in 

insertion 

locations 

date (hour, minute, second) 

101, 102 need not be set) HM.S 

121 

BOLAT 

Park orbit 

insertion 

latitude 

degrees 

122 

BOLON 

Park orbit 

insertion 

longitude 

degrees 

123 

BOVAZ 

Park orbit 

insertion 

azimuth 

degrees 

124 

BOALT 

Park orbit 

insertion 

altitude 

kilometers 

323 

324 

XIN 

Initial launch parameter values 

radians 

seconds 

325 





seconds 

326 





radians 

327 





radians 

328 





km/ sec 

343 

XLIM 

Launch control limit 

levels (see XNI) radians 

344 





seconds 

345 





seconds 

346 





radians 

347 





radians 

348 





km/ sec 
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Table 3 (continued) 


CONSTRAINT SPECIFICATION 


Location Name Meaning Units 

351 PNI 

352 

353 

354 

355 

356 

357 


PNI 

Name 

Constraint 

Desired Value 
(PSID) Location 

Tolerance (OK) 
Location 

Units f 

1 . 

X 

Cartesian state 

361 

381 

kilometers 

2. 

Y 


362 

382 

kilometers 

3'. 

Z 


363 

383 

kilometers 

4. 

VX 


364 

384 

km/ sec 

5. 

VX 


365 

385 

km/ sec 

6. 

vz 


366 

386 

km/ sec 

7. 

R 

Spherical state 

367 

387 

kilometers 

8. 

Lat 


368 

388 

degrees 

9. 

Lon 


369 

389 

degrees 

10. 

V 


370 

390 

km/ sec 

11. 

PTH 


371 

39i 

degrees 

12. 

AZM 


372 

392 

degrees 


T-7 


Constraint selector key (a number between none 
1. and 20. inclusive selects each con- 
straint, 0. ends the list - see folded 
sub -table below) 
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Table 3 (continued) 


PNI 

Name 

Constraint 

Desired Value 
(PSID) Location 

Tolerance (OK) 
Location 

Units 

-13. 

BT 

Miss vector 
components 

373 

393 

kilometers 

14. 

BR 


374 

394 

kilometers 

15. 

HEV 

Hyperbolic excess 
velocity 

375 

395 

km/ sec 

16. 

RCA 

Radius of closest 
approach 

376 

396 

kilometers 

17. 

INC 

Inclination 

377 

397 

degrees 

18. 

LAN 

Longitude of the 
asc. node 

378 

398 

degrees 

19. 

APF 

Argument of peri- 
focus 

379 

399 

degrees 

20. 

TFL 

Flight time 

380 

400 

DH.MS, sec 


Location Name Meaning Units 


419 COORD Constraint coordinate selector key none 

(1. EE 50 

2. target's orbital system, R, RxVxR.RxV 

3. target's body-fixed 

4. ecliptic and equinox) 


INTEGRATED TRAJECTORIES (HOW « 1.) 


6 

CDEQ 

Output interval (set it large) 

13 

TSTEP 

Factor for computing heliocentric 
integration step size (try .005) 

115 

EXTRA 

Extra de-bug output key (use 0.) 


seconds 

none 

none 


116 XTRREF Extra output key for rectification none 

and patch conditions (0. suppress, 

1 . include) 
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Location 


418 


Table 3 (continued) 

Name Meaning 

SCANNING OPTION 

SCAN Scanning option key (0. for normal 
search, n. for stepping the n-th 
selected control (see TYPE, XNI) TRIES 
times in steps of limit level (see XLIM) 


T-9 


Units 


none 
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APPENDIX A 

MATHEMATICAL DESCRIPTION OF THE START-UP CAPABILITY 



Figure A-l Heliocentric Transfer Geometry 

The fundamental equation solving Lambert's Problem for the heliocentric 
conic joining two (massless) planets is (A.l). 

r r (1 - cos ty ) 

p = — ^ (A.l) 

(r^ - r 2 cos t|r 4- sini|i tan ) 

In the equationj p is the orbit's semi-latus rectum, r^ and r^ the 
heliocentric radii of the launch and target planets, ijf is the transfer 
angle and y ^ is flight path angle at launch. The flight path angle, 
y is varied until p results in the correct transfer time through 
Kepler's Equation. The hyperbolic excess velocity, S, at launch is 
computed by A. 2. 

S = V x - V L (A. 2) 

where is the heliocentric transfer orbit's velocity at launch and 

A-l 
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V L is the launch planet's velocity. Launch time (if the launch body is 
Earth) is the time when S is contained in the Earth-fixed parking orbit 
plane, i.e., the time when (by Earth's rotation), 


H. S = 0, ‘ 


(A. 3) 


where H is normal to the parking orbit. The true anomaly, 0 , of the 
hyperbolic excess velocity vector, S, is given by 


cos 0 g 



where e is eccentricity, 

r 

e = i E. 

a 


(A. 4) 


(A. 5) 


where is parking orbit radius and where a is semi -major axis of the 
departure hyperbola. 


a s-** (A. 6) 

8 | 

r patch 

Knowing 0g, then, the angle and time in parking orbit are easily deter- 
mined. The radius vector, R, relative to the departure hyperbola is 
next calculated and used in the determination of the velocity, V, re- 
quired at injection to attain S. 


V 



4a 

r(l+R-S) 


1 



4a 

r(l+R*S) 



(A. 7) 


In A. 7, is - 
injection impulse, 
orbit velocity, V 


JA 

a 


, r is j R j and R and S are unit R and S. The 
i V, is the difference between V and the parking 


& V 



(A. 8) 
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The launch control parameter set is completed by computation of the 
azimuth, elevation and magnitude of AV. 

DVAZ _ sSiL_\ 

AV * (HxR) / 

DVEL = tan" 1 / A X. ; , R _ — 

\ / r? * 2 

\/( AV*R1 + (AV-H) 

DELV = [ AV | (A. 11) 

The launch control parameters may be used to re-compute cartesian compon- 
ents of R and V through subroutine START. 


(A. 9) 


(A. 10) 


The matched -conic iteration scheme uses R and V to compute the state, 

R and V , and time, t, , at the sphere of influence. The heliocentric 

ic 

radius vector, , is the sum of R and the launch body's heliocentric 
position at patch time, . 


Rj* = R* + Rj^ (A. 12) 

At the other end of the heliocentric trajectory, the desired values of 
B*T and B*R are used in computing radius of closest approach and eccen- 
tricity of the arrival hyperbola. These permit computation of the time, 

* 

t^ , of patch to the target's sphere of influence and the target body's 
heliocentric position at that time, R_. The hyperbolic excess velocity 

•** A a 

at arrival, S, enables computation of unit vectors, T and R, in the miss 

plane. The desired miss-vector, B, and the target -centered position 
* 

vector, R , at target-patch are then computed. 
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B = (B«T) T + ( B*R) • R (A. 13) 

R = - sin (0+a )B + cos ( 0 + 0 ' )s (A. 14) 

In A. 14, 9 is the negative true anomaly at patch on the arrival hyperbola 
and a is the half-angle between the asymptotes. The heliocentric position 
at patch, R^ , is computed by A. 15. 

R 2 * * R* + R t (A. 15) 

^ ^ 

By using R^ , R 2 , t^ and t instead of R^, R 2 and the original departure 
and arrival dates, the iteration on equation A.l finds the heliocentric 
trajectory between spheres of influence (patch points). This trajectory 
provides improved estimates of departure and arrival hyperbolic excess 
velocities which may then be used to calculate launch control parameters 
and perhaps initiate another conic-matching iteration. 
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APPENDIX B 

MATHEMATICAL DESCRIPTION OF THE SEARCH CAPABILITY 


The control variables for the search capability are selected from one 
of the three sets shown below. 


AVAILABLE CONTROLS 


Cartesian 

Spherical 

Launch 

Parameters 

X 

R 

DLAZ 

Y 

LAT 

DTL 

Z 

LON 

PRKT 

vx 

V 

DVEL 

VY 

PTH 

DVAZ 

VZ 

AZ 

DELV 


The initial conditions of the trajectory are completely specified by any 
one of these sets appropriately transformed into cartesian EE50 compon- 
ents. Let S(t) represent the cartesian EE50 state vector at time t, 
let X i represent the control set for the i-th trial and let Y (a constant 
vector) complete the set required to compute S(t Q ) . That is, 

S i (t Q ) = F(X r Y) (B.l) 

where F represents the transformation which maps X^^^jY into S at t Q . 

The transformation is performed by subroutine CONVX for the cartesian 
or spherical sets and by START for the launch control set (see Reference 
2 for details) . 

The state at t is derived from S(t Q ) by patched-conic trajectory calcu- 
lation or by numerical integration of the perturbed equations of motion. 
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Let this trajectory generation function be denoted by G. 

S(t) - G(S(t o )) (B.2) 

The final time, t, may represent the occurrence of a significant trajectory 
event such as closest approach to the target as well as a specified time 
after t Q . See Appendix B of Reference 1 and Appendix C of this report 
for details relative to equation B.2. 

The constraint vector, (f , is a function of S(t). 

V = P(S(t)) (B.3) 

The components of are user-selected from an. available set of 20 
constraint functions. Subroutine ENDCON computes the relationship, P, 
indicated by B.3. The constraint functions are related to the controls, 

X, through the transformation or ,, plant n indicated by equations B.l, B.2 
and B.3. The gradient, H, which represents (linearly) the sensitivity of 
if to changes in X is computed by the method of finite differences. 

H J „ _3l_ « Ks+axb - Ifia (b. 4 ) 

*X J AX j 

Equation B.4 symbolically represents the method of computing the j-th 
column of H which is the sensitivity of ijr to the j-th control variable, 

X" 3 . X J is incremented by X" 3 to form a new control vector, X+ X J . 

This control vector is used, through B.l, B.2 and B.3, to compute a new 
t (X+ AX' 3 ). The sensitivity, H' 3 is then computed as in B.4. When H 
has been completed by computation of each of its columns, it is used 
as shown in B.5 to compute the control set for the (i+l)-th trial- 

T T " 1 

X i+1 = X i + H (HH > ( X i>^ <».5) 

In B.5, ’jf D is the desired constraint vector as specified by the user. 


B-2 


PHILCO 


SPACE G. RE-ENTRY SYSTEMS DIVISION 
Philco-Ford Corporation 



TR-DA2154 


The available constraints are shown below. 

AVAILABLE CONSTRAINTS 


Cartesian Spherical Other 


1 

X 

7 

r 

B 

B*T 

19 

2 

Y 

8 

lat 

B 

B‘R 

20 

3 

Z 

9 

Ion 

15 

v <» 


4 

VX 


V 

16 

r 

P 


5 

VY 

11 

Y 

17 

i 


6 

VZ 

12 

az 

18 

n 



Most of the available constraints are defined relative to some specific 
coordinate frame. The coordinatization is implied in B.3 and must be 
uniform for all elements of i|/ . Available coordinates include Earth's 
mean equator and equinox of 1950.0, target body orbital, target-fixed 
and ecliptic-equinox. 
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APPENDIX C 

NUMERICAL INTEGRATION 

This appendix describes the mathematical theory of numerical integration 
used in the Mark IV Error Propagation Program. It is implemented in sub- 
routine DEQS. This appendix is taken directly from Reference 4 and still 
contains equation sequence numbers from that document. 


We now consider the integration of the equations of motion and variational equations 
for a given set of parameters, U, and given initial conditions R(t^), V^). All 
the equations may be considered as the vector equation 

*X = f(X, X) (D. 4-1) 

In any numerical integration process, we approximate the integralX(t) at a sequence 
of points, t., on the integration interval, (t Q , t^) obtaining the X (t.) from some 
approximation of the Taylor* s series 

X (t i + 1 X(t.) +hX(t.) + |h 2 X(t.) +|h 3 x‘(t.) +... 

h = fc i + 1 _t i (D. 4-2) 
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At any t, the second derivative may be determined from the differential equation, 
and the higher order derivatives must be developed implicitly from the known 
derivative at neighboring points. The various methods differ in the way in which 
the series (D. 4-2) is approximated. 

The ODP uses Adams* method, which approximates the series using the values 
of f (X,X) computed at the previous integration points, t. t. , etc. , for 

1 “ i. 1 *“ Ct 

long-term integration, and a generalized Kutta method for short-term integration 
and for starting the Adams’ integration. The Kutta method uses values of f (X,X) 
at suitably chosen points on the internal (L, L + ^). The two methods are des- 
cribed below. 

D. 4, 1 Adams 1 Method 

We assume that the quantities 


X. 

1 

- X(t.) 


X. 

1 

- X(t.) 

(D. 4-3) 

f. 

1 

= f(X r X.) 



have been determined at the sequence of equally spaced points 


t 

n - m 


t - mh 
n 


m = 0, 1 N 


(D. 4-4) 
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We write the Taylor's series 

f <t n + sh) - f n + sh f(t n ) + | s 2 h 2 f (y + . . . (D. 4-5) 

N 

truncating after terms in-(sh) . The coefficients of the resulting Nth degree 

» 

polynomial is s may be determined to satisfy the N + 1 conditions. 

f (t - mh) - f (D. 4-6) 

' n ' n - m ' ' 


The polynomial is usually written in terms of the backward differences 




1 



= Vf - 
n 


Vf. 


n 


1 




f 

n 


1 


and hence 


f (o) 

n + s 


N 


£ 


a k (s)? k f n 


a =1 
o 

a fc = s (s + 1) . . . (s + k - 1), k > 1 


(D. 4-7) 


(D. 4-8) 
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The error in approximation on the interval (L, t . + is 


<°> = a (s) h N + 1 f (N+1) (S) 


f v sh >- f r; s = a N + i 


(D. 4-9) 


t -N s; § < t + sh 
n — n 


D. 4. 1. 1 Integration Formulas. If we substitute the polynomial (D. 4-8) into 
the integral relationships 


t +sh 
n 

X n + S “ X a *J f(x (t), X(t)) dt- 

t 

n 


t +sh t 
n 1 


X n + s = K + f f f (X (t), X (t>) dt dt (D. 4-10) 


t t 
n n 


we obtain the approximations 


. (o) 
n + s 


N 


x „ + h £ v s > vkf n 

k = 0 


(o) 

r 

““n + s 


X + shX + h 
n n 


N 


k' = 0 


E k (s)V f n (D.4-11) 
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where 


A k (s) = f \ (t)dfc 


B.(s) 


s 

-! 


■'j c \ °/ - j (t) dt 

o 


(D. 4-12) 


with errors 


X - X =A / o\ u N + 2 ^ (N + 1) .p . 

n + s n + s A N + l< s > h 1 <V 

X. -X <?> = (s) h N + 3 f ( N + ^ (Tlo) (D. 4-13) 

n+s n+s N + 1 w v ' v ' 


t XT < § , %> < t + sh 
n - N- o - n 


since + ^ (s), + (s) do not change sign on (0, 1). These formulas resulting 

from extrapolation are termed open . An alternative form of the polynomial 


N 


f = ]P a, (s— 1) V k f ,, 
n + s a— ' ' k v ' n+1 


(D. 4-14) 


k=o 


yields the closed formulas 


N 


X £> = 
n + s 


X n + h Z C k < s > V k f 


k=o 


n+1 
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X d) 
n+s 

C k (B) 

D k (B) 


N 

= X + shX +h 2 2 D J S ) V k f 
n n .k n+1 

s 

“/ a k (t-1) dt 
o 

s 

=/ C k (t)dt 
o 


with errors 

Vs-^Ss - C N+1 (s) h N+2 ^(S,) 

X n +S - X iil - Vl <S) h N+S i< N+1 H> 


(D. 4-15) 


(D. 4-16) 


The closed formulas require knowledge of f(X^ + , , X j) for the determination of 
^n+r ^n+1’ and ^ ence ma y be used directly only in simple quadrature. For the 
integration of differential equations, they must be used in conjunction with formulas 
for the prediction of X^ , The obvious solution is to use the open formulas 

as predictors to compute estimates xj^, x|°^, and to use the closed formulas 
as correctors. Evaluating the coefficients at s = 1, 


X = x (< 2 
n+1 n+1 


+ A h N+ 2 f {N+l) ^ j 
n+1 ' ^ 


X 


( 1 ) 

n+1 


N 


+ C N+1 h ' 


N+2 
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X 


n+1 


X ( °> + B„ h N+3 f ( N+1 ) ( T1 ) 

n+1 N+1 ' o' 


N 


-(D. 4-17) 

(o) 


= xW ♦ D Nn h™ * T, D/(f n+1 - IW) 


If we assume that h is sufficiently small that h(f n+ ^ - 4+j) is negligible compared 
with h N+2 f^ N+1 ^ (§), and that / N+1> (S) varies only slowly with §, we may elimi- 
nate X +1 , X , obtaining 


h N+2 j(N+l) @ = ^ . X W ) /(A n+1 - C N+l ) (D. 4-18) 


Using the easily established relations 


s ) ^k+1^ "* ^t+i^ 


k+1' 


-.k- _ k,. k+1, 

V f n V f n+l " V f n+l 


(D. 4-19) 


v/e have 


h N+1 f< N+1 >(S) V N+ \ +1 


(D. 4-20) 


and hence our best estimate of the integrals is 


X_,,' = X ( _°\ + hAxrj-i V N+1 4+1 


n+1 n+1 


N+1 


- X<°i, + h 2 B N+1 , N+1 a 


n+1 n+1 


n+1 
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The integration coefficients are listed through k=8 in Table D-2, below. 


TABLE D-2 

ADAMS' INTEGRATION' COEFFICIENTS 


k 

A k 

B k 

c k 

D k 

0 

1 

1 

2 

1 

1 

2 

*t 

1 

1 

1 

1 

± 






2 

6 

2 

3 

o 

5 

3 

1 

1 


12 





24 

12 

24 

3 

9 

38 

1 

7 ' 


24 

360 

24 

- 360 

4 

251 

135 

19 

17 


720 

1440 

720 

1440 

r 

• 475 

863 

27 

82 


1440 

10080 


10080 

6 

19087 

9625 

863 

731 


60480 

120960 

60480 

120960 

7 

36799 

135812 

1375 

' 8563 

i 

120960 

1814400 

120960 

1814400 

8 

1070017 

515529 

33953 

17719 

3628800 

7257600 

3628800 

7257600 


D. 4. 1. 2 Interpolation . To obtain f, X, X at points other than integration points, 
we may use the polynominals (D. 4-8), (D. 4-11). Setting 


t - t + sh 
a 


_ . h k d k f(V . 
h ^ - 


d k f<y 

7E 

ds 


(D. 4-22) 
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we obtain 


N 

f (t) =.J2 F k s k/k I 

k=o 


N 


k+1 


X(t) = X n + h > F r s /(k+1) I 


N 


X(t) = X n + sh X n + h 2 2 F k s k+2 /(k+2) J 

k=o 


and for s on the interval (-1, 0), the derivatives are obtained from 


— — 


— 







— 


— - — 

F o 


1 

0 

0 

0 

0 

0 

0 

0 


f 

n 

F 1 


0 

1 

1 

2 

1 

3 

1 

4 

1 

5 

1 

6 - 

1 

7 


V f n 







11 

5 

137 

7 


V 2 f 
v n 

j 

F 2 


0 

0 

1 

1 

12 

6 

180 

10 








3 

7 

15 

29 


v 3 f 

v n 

F 3 


0 

. 

0 

0 

1 

2 

4 

8 

15 


F 4 


0 

0 

0 

0 

1 

2 

17 

6 

7 

2 


v 4 f 

v n 

F 5 


0 

0 

0 

0 

0 

1 

5 

2 

25 

6 


v 5 f 

n 

F 6 


0 

0 

0 

0 

0 

0 

1 

3 


<1 

CD 

s" 1 

F 7 


0 

0 

0 

0 

0 

0 

0 

1 


7 

v f 

v n 











— 


_ — 


til 

where all differences after the N are to be set' zero. 
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D. 4. 1. 3 Change of Interval Size. For a set of differences V * n for the spacing h, 
we may compute an equivalent set for any spacing sh, so that the interpolation 
polynominala for the two sets are identical In t. Two particular changes may be 
made rather simply, for s = l/2 and s = 2, and these changes provide all the 
spacing flexibility required. 

Using (D. 4-22), we have for s = 1/2, 
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and for S - 2, 


— — 


— 







f 

n 


1 

0 

0 

0 

0 

0 

0 

V^n 



2 

-1 

0 

‘ 0 

0 

0 

v 2 f 

n 




4 

-4 

1 

0 

0 

v 3 f 

n 





8 

-12 

6 

-1 

A 






16 

-32 

24 

A 



4 




32 

-80 

v 6 f 

n 








64 

V 7 f 
v n 

L— — J 


_ 








0 


f 



n 

0 


v f „ 

0 


v 2 f 

v n 

0 


v 3 f 

n 

-8 


V 4 f 
v n 

80 


V 5 f 
v n 

192 


A 

128 


/'»_ 


(D. 4-26) 


D. 4. 1/4 Ordinate Formulas. . The use of difference formulas has some computa- 
tional disadvantages. At each integration point, a complete set of differences must 
be computed, and the old set must be retained until the integration accuracy is 
verified. More efficient computation results from direct use of the computed 
ordinates. The corresponding formulas may be obtained from the relations 


k 

k f . ( -1)*" m! { 

v n / j k! (k-m)! n-m 
m=0 

/ 


(D. 4-27) 


The various coefficients depend upon N as well as on k. For N = 5, the integration 
formulas are: 


X<°> 

n+4 


X n + 10080 


29939f - 55461f , + 69874f „ 

n n-1 n-2 


- 51086f 0 + 20139f . - 3325f _ 

n-3 n-4 n-5 


C-ll 
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X (°) 

n+l 


X n + 15580 


10852f - 15487f , + 18752f _ 

n n-1 n-2 


13474f „ + 5260f . - 863f c 

n-o n-4 n-5 


V - X <°) 4 - 19a87h T 7 6 f 

n+l ~ n+l 60480 v n+l 


X - t(°) + 9625n 6. 

n+l “ n+l 120960 v n+l 


v 6 f 


= f (0 ! - 6f + 15f , - 2 Of . +15f „ - 6f „ + f 


n+l n+l 


n 


n-1 


n-2 


n-3 


n-4 n-5 


(D. 4-28) 
(Contd. ) 


The interpolation formulas are 


F 

o 


60 0 

0 

0 

0 

0 


f 

n 

F 1 


137 -300 ' 

300 

-200 

75 

-12 


f n-l 

F 2 


225 -770 

1070 

-780 

305 

-50 


f n-2 

F 3 

1 

60 

255 -1065 

1770 

-1470 

615 

-105 


£ n-3 

F 4 


180 -840 

1560 

-1440 

660 

-120 


f n-4 

_ F 5_ 


60 -300 

600 

-600 

300 

-60 


f n-5 


(D. 4-29) 


D. 4. 2 Generalized Kutta Method 


The various methods called Kutta or Kunge-Kutta methods are based on a process 
suggested by Runge (Reference 7) and developed for first order equations by Kutta 
(Reference 8), Applied to second order equations, the method requires evaluation 
of the derivative f at the sequence of points : 
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t 

n 


t + 
o 


C h 
n 


n-1 


X 

n 




.f. 

ni i 


X 

n 


X + h 


b X 
\ n o 


n-1 

i=0 


d .{. 

ni i 


C 

n 


i=n 


(D. 4-30) 


where t is the initial point on the integration interval. The process (D. 4-30) 
is repeated through N substitutions (n = 0,1,. . . , N-1), and X(t -*-h),‘ X(t Q +h) are 
then approximated by the N + 1 st values in the sequence. Appropriate values of 
the coefficients are determined by matching as many as possible of the leading 
terms of the series (D. 4-2) with those of the series obtained by substituting the 
Taylors series for f(X,X) into the-sequence (D. 4-30). 


Several methods have been developed for integrating (D. 4-1) and for using two 
adjacent intervals for computing truncation error, interpolating between interval 
end-points, etc. Miachin (Reference 9) treated the special case X = f(X), obtaining 
accuracy through terms in h and an expression for the truncation error in 
X(t +h) using derivatives computed on the two intervals (t , t Q +h) and (t Q +h, t Q +2h). 
In two unpublished communications, T. W. Hinton treated the case X = f(X,X), 
obtaining accuracy through h for the case df/dX = 0 and through h for the general 
case. Hinton also gave an expression for 'the truncation error in X(t Q +h) and 
equations for interpolating on the interval pair (t , t +h) and (t Q +h, t Q +2h). 


* Lockheed Missiles and Space Co. , Sunnjrvale, California, August 1963. 
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D.4. 2.1 Integration Formulas. The coefficients given by Hinton are: 


C x = b x = 3/10 
C 2 = b 2 = 3/4 


C 3 = b 3 “ C 4 = b 4 " 1 


dj = 9/200 


d = 9/32, 


d 3 - d 4 -l/2 


20 ' - 21/32 

d 20'° 

21 = 45/32 

d 21 - 9/52 

30 = 83/27 

d 3Q - 10/27 

31 = -280/81 

d 31 = 7/162 

32 = 112/81 

d 32 = 14//81 

40= 5/54 

d 4Q = 5/54 

41 = 250/567 

d 41 = 25/81 

42 = 32/81 

d 42 = 8/81 

43 ' ^ 14 

d 43 = ° 


C 43 = 1/14 d 43 = 0 (D. 4-J 

If we denote by f. and f. 0 the values calculated for f (X.,X.) on the Intervals 
1)1 1 > ^ 11 

(t Q , t +h) and {t Q +h, t +2h), respectively, the truncation error in X(t Q +2h), 
assuming no error in X (t Q +h) is 

i2 r 


T ~ 34020 81f : 


81f 3,2 + 112f 2,2 - 550f l,2- 2478f C 


+ 1134 ^ + 3248 ^ 


2450f l,l + 903f 0>1 


(D. 4-32) 
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5 6 

This term is of order h and the error in the approximation is of order h for 

the general case. 


D, 4. 2. 2 Interpolation. Linear combinations of the f may be used to interpolate 
for f,X, X on the interval (t o ,t o +2h). We again set 



(D. 4-33) 


obtaining the interpolation formulas (D. 4-23), The are given by 


r-.i 


486 

1344 

-1200 

1638 

-486 

-1344 

1200 -504 



1620 

2912 

-8900 

8904 

-2106 

-5600 

5900 -2730 ' 

— F 
2 2 

1 

1458 

2016 

-9900 

9828 

-1944 

-4704 

6900 -3654 

1134 

— F 
6 3 


324 

448 

-2200 

1428 

-324 

-448 

2200 -1428 


3,2 

f 2,2 


1,2 


0,2 


3,1 


2,1 


1,1 


0,1 


(D. 4-34) 


3 4 5 

The expression yields accuracy through h ,h .,h for f, X, X respectively for the 
general case f(X,X). 
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D. 4. 2. 3 Conversion to Adams 1 Ordinates . As we noted earlier, some starting 
process is required to accumulate the necessary ordinates for the Adams' 
integration. The necessary ordinates are computed by the ODP by interpolation 
on a single interval pair integrated in the Kutta mode. To avoid extrapolation 
beyond the interval pair in computing V f n » the highest significant difference 
obtainable, we take S = 1/2 and set V 1 ^ - 0 tor all j ^ 4. We have 



896 

-800 

1092 

-324 

-896 

800 

672 

-1500 

1449 

-405 

-1120 

1000 

44B 

-2200 

2562 

-486 

-1344 

1200 

224 

-1100 

714 

-162 

-224 

1100 



(D. 4-35) 
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APPENDIX D 
INPUT SUMMARY 


I CARD FORMATS 


D.1.1 General Formats 


5 



EXEC 
1 = 1 

PINT 

ERP 

SPOUT 

IM 


1 • 

— 2 




= 3 




= 4 

CONW 




= 5 

START-UP 




= 6 

SEARCH 




= 0 

terminate run 
and wrap up tapeb 

3 

61 

72 

HEADER 

Cards 




1 k| 



1 HEAD 












Any 

alphanumeric information 



3 


18 

61 





11 

IN 

i 

I 



6 

12 18 

72 

11 , '"1 


12 





ini 

1 l 



K - 0 or 1 


Integer cards 

11 = 1st subscript 
IN *= end subscript 


15 18 30 33 45 48 60 61 


Ej 


12 

V2 

a 

V3 

a 

V4 | Unit 1 


(4(I3,E15.8)A6) 


OVERIAY cards 

Ij_ = subscript for 

V. = value, 
x 

Unit = FT 

FT/SEC 


NM 

NM/SEC 


SQUARE 
(or blank) . 
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D. 1.2 


PHILCO 


Special ERP Formats 

Overlay card with II - NCH 

Vl “ TSTART 

V2 = TSTOP 

V3 = OINTV 


NCR Control. times card 

NCH SO, process measurements from 
TSTART to TSTOP. 

NCH >0, read processing options or 
changes . 

NCH =0, no changes to read. 

NCH <0, no measurements to be 
processed . 

NCH=111, stop. 


1 7 9 

11 13 

15 

|N[ NAME 

h 


H 

*4 

1 I 


3 [ 


55 57 

. Change card 

□ 

i 24 

i 25 

U 

j n = type of change 


1, stations 

2, beacons 

3, onboard 

4, equation of motion 


i, = error source 

treatment code. 

(See ERP section, tables 
1, 2, 4, 6, 8, 10). 


Overlay Card 


Data Cards 

(Station, beacon, etc.) 

(See ERP section, tables 
3, 5, 7, 9). 
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II INPUT ORDER AND DESCRIPTIONS 
D.2.1 PINT 

Input for a PINT case consists of: 

1. EXEC card (1 in column 5) 

2. C-array data, ending with blank. This includes the header 
for the vehicle ephemeris tape both when reading and writing 
this tape (K = 0 format) . 

3. IN, X-array data ending with blank. 

Both blanks must be included, whether or not any data is being read 
in. 


INPCOM (C-ARRAY) PARAMETERS RELEVANT TO PRECISION INTEGRATION 


Address 

Name 

Dimension 

Description 

C(8) 

ASTU 

l' 

Conversion factor, kilometers/a .u. 

C<21) 

UM 

10 

3 2 

Gravitational constants (km /sec ). 

C(31) 

RPL 

10 

Planetary semi-major axes or radii 
(km), Earth's radius value used in 
gravity and park orbit calculations. 

C(51) 

RPAT 

10 

Sphere of influence radii (km) or 
distance at which transfer is made 
to or from each body. 

C(71) 

WE 

1 

Earth's sidereal rotation rate 
(rad/sec) . 
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C(585) EMP 


C(609) ZH 


C(613) XH 


C(623) ABC 

C(626) DRAGC 

C<628) TMAG 

C(629) SPK 

C(630) DUTD 


24 


4 


10 


3 


2 


1 


1 


1 


Indicators for equation of motion 
parameters to consider for computing 
sensitivities. These indicators 
are set 1. in INPCOM block data 
A zero in EMP(l) asks for the 
sensitivity of state to EMP(I), 


Zonal harmonic coefficients, of 
the Earth's gravitational field 
(dimensionless) . 


The TH array contains 10 tes serai 
harmonics ordered J 21 , X 21 » J 2 2* X 22* 

J 31* X 31* J 32* J 33’ X 33* 
til 

J is the m tes serai harmonic coef- 

nm th 

ficient of the n order, (no units); 

and X is the geographic longitude 
nm 

corresponding to J . 

nm 

Mass- normalized moments of inertia of 

the moon used for computing the triaxial 

2 

gravity perturbation (km ). 


Atmospheric drag coefficients of the 
model -c 1 e“ c 2 h | V I V. Both have 
units of (km ). 


Magnitude of thrusting acceleration 
(km/sec^) , 

Solar radiation pressure coefficient 

(knfVsec^). 

Discrepancy between ephemeris time and 
universal time (days). 
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C(672) XCASE 1 Case counter or trajectories generated 

and written on the vehicle ephemeris 
tape by a PINT run. Set initially zero 
by INPCOM BLOCK DATA, ICASE is incre- 
mented automatically and written on tape, 
so that subsequent tape-uBing routines 
may select a desired trajectory by this 
identification number. 

C(675) HEAD 12 Alphanumeric header written on the vehi- 

cle ephemeris tape, and used by reading 
routines to identify which tape to read. 


Address 

IN(1) 


IN(2) 
IN (3) 


FIXED-POINT WCCM (IN-ARRAY) PARAMETERS FOR THE 
PRECISION INTEGRATION OPTION 


Name 

IB1 


IB2 

KREF 


Dimension Description 

1 Launch body number or initial central 

body (1 Earth, 2 Moon, 3 Sun, 4 Venus, 

5 Mars, 6 Saturn, 7 Jupiter). 

1 ' Intended target body number. 

1 KREF < 0 signals no refine, but compute 

the state transition matrix. 


KREF = 0 signals no refine or state 
transition matrix, but trajectory 
calculation and maybe tape-write. 


KREF = 1 signals refinement option 
but do not save the solution. 

KREF = 2 signals refine and save. 
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IN (4) 

NRIPBO 


Tape-write or sensitivity computation 
indicator (0) for neither^ (1) for 
writing a precision trajectory on tape, 
(-K) for integrating variational equa- 
tions from case #K on tape for motion 
parameter sensitivities. 

IN (10) 

NDEQ 

9 

Integration package indicators, (see 
DEQ writeup for definitions--need not 
be changed for normal operation.) 

IN(19) 

IBC 

10 

Bodies to be considered as perturbing 
force centers , (1) consider^ 0) omit. 

IN (29) 

NH 

5 

Tesseral harmonic gravity indicators 
NH(1) is the order of the highest 
zonal harmonic to be used (for J use 5), 

NH(2-5) is the highest degree longitudinal 
harmonic to be included for each order 
(NH(3) *» 2 means use and J^)- 


D-6 


PH ILCO 


SPACE S. RE-ENTRY SYSTEMS DIVISION 
Philco-Ford Corporation 



TR-DA2 154 


Address 

X(l) 


X(9) 

X(10) 

X(ll) 


FLOATING-POINT WCOM (X- ARRAY) PARAMETERS FOR THE 
PRECISION INTEGRATION OPTION 


Name 

CDEQ 


RECT 


RN00B 


TSTEP 


Dimension 

8 


1 


I 

10 


Description 
Integration parameters 

CDEQ(l) ~ current step size computed 
by program. 

CDEQ(2) *= first (next) time at which 
output is desired (sec), 
reset each case. 

CDEQ (3) = initial step size, set by 
program. 

CDEQ (4) = doubling limit, set by 
program. 

CDEQ (5) = halving limit, set by program. 

CDEQ(6) ° output time interval (sec) 
starting at CDEQ(2). 

CDEQ(7) = minimum divisor for relative 
error. 

CDEQ(8) c upper bound allowed on 

truncation error, run time 
increasing as CDEQ (8) is 
decreased. 

Conic rectification tolerances. Reference 
conic section for Encke's method is 
changed when ratio of position deviation 
to radius from central body exceeds RECT 1 . 

Radius from the Earth at which oblate- 
ness ceases to be considered (km). 

Fraction of a radian to use in computing 

the initial integration step size 

( ""N n 

— . rounded to 2 for 

v y 

each central body). 
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X(21) ERRTBlL 1 

X(22) XTRTAP 1 

X(23) DTPREC 1 

X(24) DTNUT 1 


Relative error tolerance for inter 
polating when writing vehicle 
ephemeris tape. 

Key for extra output when writing 
ephemeris: (0.) no extra output* 

(1.) extra output. 

Time interval over which Earth's 
precession matrix may be assumed 
constant (sec) , 

Time interval over which Earth's 
nutation matrix may be assumed 
constant (sec). 
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NORMAL TRAJECTORY INPUT DATA 


Address 

Name 

Dimension 

Description 

X(101) 

YW 

1 

Year, month and day of injection 
(YM.D, e.g. 6311.16 is Nov 16, 1963.) 

X(102) 

YF 

1 , 

Hour, minute and second of injection 
(HM.S . e.g. 1408.1409 is 14 h 8 m 14.09®). 

X(103) 

B0DN0 

1 

Body number to which the injection 
state is relative (1. Earth, 2. Moon, 
etc) . 

X(104) 

TYPIN 

1 

Input state coordinate type in the 


2 

form A. 10 + B.10 + C where 

A, B and C are interpreted as follows: 

A = 0, Cartesian; A = 1, Spherical, 

A = 2, Orbital Elements 1 , B 53 0, Equator, 

B * 1, Ecliptic; B = 2, Body-fixed; 

C =.0, Epoch of 1950.0; C = 1, Epoch 
of Date (101. means Spherical, Equator 
and Equinox of Date) . 

X(105) RZ 3 Cartesian position vector (km) if 

TYPIN < 100. Radius (km), latitude (deg) 
longitude (deg) if 100. < TYPIN < 200. 
Orbital elements if TYPIN > 200: Semi- 
major axis (km), eccentricity, true ano- 
maly (deg). 
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X(108) 

VZ 

3 

Cartesian velocity vector (km/sec) if 
TYPIN < 100. Spherical velocity if 100 
£ TYPIN < 200: Speed (km/sec), path 

angle (deg), azimuth (deg). Orbital 
elements if TYPIN s 200: Longitude of 

the node (deg), inclination (deg\ argu- 
ment of per laps is (deg). 

X(lll) 

TSECI 

1 

Trajectory starting time from YW and 

2 -2 
YF. (Days. 10 + Hours + Minutes. 10 

-4 t 

+ Seconds *10 ). 

X(112) 

TSTP 

1 

Trajectory stopping time from YW and 
2 ,,,0 
YF. (Days. 10 + Hours + Minutes. 10 

+ Seconds. 10 ^) , 

X(I13) 

RSTP 

1 

Trajectory stopping radius from IB2, 
the target body (km) . 

X(114) 

OYTP 

X 

• 1 

Trajectory output indicator: 

(0.) minimum output, (1.) Equator and 
Equinox of 1950.0,(-1.) Equator and 
Equinox of date. 

X(115) 

XTRA 

1 

Extra output key for stopping 
functions in FSUB and DERIV : 

(0.) omit, (1.) print extra output. 

X(116) 

XTRREF 

1 

Extra output key for REFINE infor- 
mation:^.) none»(l.) print extra 
information. 
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“V 

REFINE INPUT DATA 
Parking Orbit Data 


Location 

Name 

Description 

X(120) 

STARTK 

(1.) Asks for "Earth-fixed" park orbit, 
(2.) Asks for "Inertial" park orbit. 

X(121) 

DLAT 

Insertion latitude (+ 90 deg) measured positive 
northward from the equator of insertion epoch. 

X(122) 

DLON 

Insertion longitude (+ 180 deg) measured positive 
eastward from the X-axis or Greenwich Meridian, 

X(123) 

DVAZ 

Insertion velocity azimuth (+ 180 deg) measured 
positive clockwise from local north at insertion. 

X<124) 

PALT 

Altitude at insertion and of the circular park- 
ing orbit (km). 

X(125) 

ORBN 

Number of whole parking orbits after insertion 
and before injection (floating-point integer). 

X(128) 

TWINS 

Year, month and day of park orbit insertion 
(format YM.D, two digits each), for example 
6608.03 is Aug. 3, 1966). 

X(129) 

YFINS 

Hour, minute and second of parking orbit inser- 
tion (format HM.S; for example 1348.0533 is 
13 h 48* 5.33 S TJMT). 




Control Parameters 

X(130) 

TL 

Insertion time increment from nominal insertion 



date (sec). 

X(131) 

PARK! 

Time in parking orbit before injection (sec). 

X(132) 

AZIM 

Azimuth of the injection velocity impulse (rad) 


relative to the parking orbit plane at injection 
measured counter-clockwise about the injection 
radius vector. 
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Location 

Name 

Description 

X(133) 

PTH 

Elevation of the injection velocity impulse 
(rad) measured' up from the horizontal at 
injection. 

X(134) 

DELTAV 

Magnitude of the injection velocity 
impulse (km/sec). 



Control Limit Levels 

X(140) 

TLMAX 

Insertion time, maximum increment (sec). 

X(141) 

PMAX 

Park time, maximum increment (sec). 

X(142) 

AZMAX 

Injection impulse azimuth, maximum increment (rad) . 

X(143) 

ELMAX 

Injection impulse elevation, maximum increment (rad) 

X(144) 

ENMAX 

Post-injection energy (vis-viva) , maximum increment 
(km^/sec^) . 

X(145) 

TRY 

Number of trials (iterations) allowed in REFINE 
if convergence is not obtained earlier. 

X(146) 

SX 

Fraction of control limit levels to be used for 
computing numerical partials. 

X(147) 

SCALE 

Distance to which end point variational positions 
are to be scaled before using the variational 
state to compute constraint partials (km). 



Desired Constraints 

X(150) 

FLTXME 

Indicator: (0.) neither flight time nor target- 


' 

relative energy are constrained, (1.) constrain 
flight time,(-l.) constrain arrival energy. 

X(151) 

TFL 

Desired time of flight to closest approach 
or impact from park orbit insertion (Days Hours. 
Minutes Seconds) . 

X(152) 

C3D 

Desired energy (Vis-Viva) relative to the 



2 2 

target body at arrival (km / sec ). 
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Location 

Name 

Description 

X(153) 

TLA.T 

Target vector latitude measured positive toward 
the z-axis of the target coordinate system 



(+ 90 deg). 

X(154) 

TLON 

Target vector longitude measured positive toward 
the y-axis from the x-axis of the target coordi- 
nate system (+ 180 deg) . 

X(155) 

PRD 

Periapsis radius desired at radius of closest 
approach (km) . 

X(156) 

PACTIM 

Indicator: (0.) no impact of the target body, 
(1.) impact the target, ignoring PRD above. 

X(157) 

HOWB 

Indicator: constrain 

-1., only flight time or energy (see FLTIME) 
0., TIAT , TLON, PRD, see FLTIME 

1., PRD, see FLTIME 

2., BDTD, BDRD, see FLTIME 

X(158) 

BDTD 

Desired B.T component if H0WB=2. (km) 

X(159) - 

BDRD 

Desired B*R component if HOWB^. (km) 


Earth Return Constraints 


X(160) 

DLATE 

Latitude of the desired touchdown point 
(± 90 deg) . 

X(161) 

DLONE 

Geographic longitude of the desired touchdown 
point (+ 180 deg) . 

X(162) 

DVAZE 

Desired inertial velocity azimuth at the touch- 
down point, measured positive clockwise from 
north (+180 deg) . 

X(163) 

PERH 

Altitude of virtual perigee desired (km). 

X(164) 

AVETIM 

Reciprocal average angular rate from virtual 
perigee to touchdown (sec/rad). 

X(165) 

RMIN 

Indicator to choose between two possible range 
angles from entry to touchdown (+min, 0 max) . 
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Location 

Name 

Description 

Constraint Satisfaction Criteria 

X(170) 

VOK 

Injection velocity impulse being minimized is 
close enough if within VOK of minimum (km/sec) 

X(171) 

BDTOK 

B*T computed from target input constraints 
and that of the test trajectory must be closer 
than BDTOK (km) . 

X(172) 

BDR0K 

B*R computed must be closer than BDROK (km) 
to B*R of the test trajectory for convergence. 

X(173) 

FLT0K 

Flight time satisfaction criterion (sec). 

X(174) 

EN0K 

Energy (vis-viva) satisfaction criterion 
(km^/sec^) . 

X(175) 

EPSM 

Factor for testing convergence in FNDMXN. 

X(176) 

VST 

Starting improvement when minimizing velocity 
.impulse magnitude (km/sec). 

X(177) 

PARKK 

Indicator for inertial parking orbits: 


(1) minimize velocity impulse by varying park 
time, (0.) fix park time. 
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2.2 ERF 

Input for an ERP case consists of: 

1. EXEC card (2 in column 5). 

2. C, IW data, ending in a blank. This includes 2 headers, 

K - 0 for the vehicle ephemeris tape and K = 1 for the case 

heading(and special output tape if requestec^.^ C(699) should 
be set ^ 0 if a tape is not wanted.^ 

The blank following this data must appear. 

3. NCH control times card. 

If NCH = 111, the next card will be a new EXEC card for the next 
case. 

If NCH < 0, the next card will be a new control times card, unless 
TSTOP > the flight time of the run. In the latter case, the next 
card will be a new EXEC card. 

If NCH > 0 (but not = 111) the next card is a change card. 

4. Change card. (Station, beacon, etc). 

4a. Appropriate .(station, beacon, etc) data cards, ending with 
a blank. This blank is necessary whether or not there are 
any data cards . 

The combination of change and data cards may be repeated as 
often as necessary. The sequence must be terminated by a 
blank to signify the end of the change cards. 
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If NCH > 0 but no changes are desired, then the deck will 
appear : 


NCH Control times card 
Blank (end of change cards) 

NCH Control times card or EXEC card as necessary. 

If changes are included, the deck may look like: 

NCH Control time card 
\ Beacon change card 
Beacon data card 
Blank (:.nd of data) 

Blank (end of changes) 

NCH control times card or EXEC card, as necessary-. 
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Code 

00 

2 

01 

-1 


Columns 

1 

2-7 

8-9 

10-11 

12-13 

14-15 

16-17 

18-19 

20-21 

22-23 

24-25 

26-27 

28-29 

30-31 

32-33 


Table 1 Error Source Treatment Code 


Meaning 

Omit any treatment of the quantity 

Treat the quantity as having only random errors (applicable 
only to measurements, i.e, , where measurement has no bias) 

Consider the quantity as a deterministic error source (and with 
random errors, if the quantity is a measurement, i.e., when a 
measurement has a bias) 

Consider the quantity as a deterministic error source in the 
sense that the error would be determined 


Table 2 Station Changes 
Quantity 

Station change indicator (1) 

Station name 

Consider or omit station (+ 1 means consider, 0 means omit) 
Station identification number (1 to 12) 

Type of angles to measured 
Range measurement 
Range rate measurement 

Azimuth, right ascension, or /-direction cosine measurement 
Elevatio^ declination, or m-direction cosine measurement 
Azimuth, right ascension, or ^.-direction cosine rate measurement 
Elevation, declination, or ra-direction cosine rate measurement 
Latitude error in station location 
Longitude error in station location 
Altitude error in station location 
Time bias error in the station clock 
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Table 3 Station 

Information 

Index Quantity 

Input Units 

2 

Period of Observation 

(Days) (Hours) . (Min. ) 
(Sec.) or (-Sec.) 
(e.g. -.5 - .5 sec) 

3 

Station Latitude (geodetic) 

Degrees 

4 

Station Longitude 

Degrees 

5 

Station Altitude 

Meters 

6 

Artifical Horizon 

Degrees 

7 

Maximum Elevation 

Degrees 

8 

Range Error(Random) 

Meters 

9 

Range Rate Error(Random) 

Meters/ Second 

10 

Azimuth Error(Random) or Right 
Ascension Error (Random) or , t- 
direction Cosine Error (Random) 

Milliradians 

Milliradians 

Unitless 

11 

Elevation Error(Random) or 
Declination Error(Random) or 
m-d±rection Cosine Error(Randera) 

Milliradians 

Milliradians 

Unitless 

12 Azinuth(etc. ) Rate Error(Random) 

Milliradians/ Sec 

13 Elevation(etc. )Rate Error(Random) 

Milliradians/Sec 

14 Range Error (Bias) 

Meters 

15 Range Rate Error (Bias) 

Meters/Sec 

16 Azimuth(etc)Error(Bias) 

Milliradians 

17 Elevation (etc.) Error(Bias) 

Milliradians 


18 

Azimuth (etc.) Rate Error (Bias) 

Milliradians/Second 

19 

Elevation (etc.) Rate Error (Bias) 

Milliradians /Second 

20 

Latitude Location Error (Northing) 

Meters 

21 

Longitude Location Error (Easting) 

Meters 

22 

Altitude Location Error (Down) 

Meters 

23 

Time Error (Bias) 

Seconds 
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Table 4 Beacon Measurement Changes 


Column 

Quantity 

1 

2 means this is a beacon change card 

2-7 

Mnemonic message 

8-9 

01 means "include", 00 means "delete" all beacons 

10-11 

00 means this is the measurement type change card 

12-13 

Number of the body on which beacons are found 
(1 - Earth, 2 - Moon, 5 - Mars) 

14-15 

Range measurement 

16-17 

Range rate measurement 

18-19 

Angle l measurement 

20-21 

Angle 2 measurement 

22-23 

Time bias error of the onboard clock 



Table 5 Beacon Measurement Change Data 

Index 

Quantity 

Input Units 

1 

Period 'of 
observations 

(Days) (Hours). (Min.) 
(Sec.) or (-Sec.) 

2 

(Spare) 


3 

Range error (random) 

meters 

4 

Range rate error(random) 

meters/second 

5 

Angle I error (random) 

miiliradlans 

6 

Angle 2 error (random) 

milliradians 

7 

Range error(bias) 

meters 

8 

Range rate error (bias) 

meters/ second 

9 

Angle 1 errors (bias) 

millirad ians 

10 

Angle 2 error (bias) 

milliradians 

11 

Time or clock' error (bias) 

seconds 
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Table 6 Individual Beacon Change Card 


Column 

1 

2-7 
8-9 
10-11 
12- 13 
14-15 
16-17 
18-19 


Quantity 

2 means this is a beacon change card 
Mnemonic for beacon identification 
01 means "include", 00 means "delete" 

Beacon number (must be 1-10) 

Number of the body on which beacon is located (same for all beacons) 
Latitude error in beacon location (Northing) 

Longitude error in beacon location (Easting) 

Altitude error in beacon location (Down) 


Table 7 Beacon Specification Data 


Index 

2 

3 

4 

5 

6 

7 

8 


Quant ity 
Latitude 
Longitude 
Altitude 

Artificial Horizon 

Latitude Location 
Uncertainty 

Longitude Location 
Uncertainty 

Altitude Uncertainty 


Input 

Degrees 

Degrees 

Meters 

Degrees 

Meters Northing 
Meters Easting 
Meters 
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Column 

l- 

2-7 

8-9 

10-11 

12-13 

14-15 

16-17 

18-19 

20-21 


22-23 

24-25 

26-27 

28-29 

30-55 


Table 8 Onboard Change Card 
Quantity 

3 means this Is an onboard change card 
Mnemonic for onboard identification 

01 means "include", 00 means "delete" onboard consideration 

Number of the body for onboard radar- type measurements 
(1-Earth, 2-Moon, 4-Venus, 5-Mars, 6-Saturn, 7-Jupiter) 

Height measurement 

Height rate- measurement 

Time bias error of the onboard clock- 

Number of angular measurements' of the first type in the cycle 
First type of angular measurement 

00 means no measurements in the stated interval 

01 means subtended angle 

02 means right ascension or longitude-type angle and 
declination or latitude-type angle 

03 means maximum llne-ef-sight change star-planet angle 

04 means minimum line-of-sight change star-planet angle 

Number of the body on which the first type measurement is made 
(1-Earth, 2-Moon, 3-Sun, 4-Venus, 5-Mars, 6-Saturn, 7-Jupiter) 

Number of angular measurements of .the second type in the 
cycle. 

Second type -of angular measurement 

Number of the body for second type angular measurement 
More specification of angular measurements 
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Table 9 Onboard Measurement Data 


Index 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 


Height error (random) 

Meters 

Height rate error (random) 

Meters/second 

Height error (bias) 

Meters 

Height rate error (bias) 

Meters/second 

Time bias uncertainty 

Seconds 

Altitude to cease radar observations 

Kilometers 

Period of radar type measurements 

(Days) (Hours) • 


(Minutes) (Seconds ) 


or (-Seconds) 

Angle 1 error k^ (random) 

Arc seconds 

(of error model = k^ + k^Zsin ^ ^)^) 


Angle 1 error k^ (random) 

Unitless 

Angle 2 error k^ (random) 

Arc seconds 

Angle 2 error (random) 

Unitless 

Subtended angle error k^ (random) 

Arc Seconds 

Subtended angle error (random) 

Unitless 

Period of angular measurements 

(Days) (Hours) • 


(Minutes ) (Seconds ) 


or (-Seconds) 

Right ascension of reference star 

Degrees 

Declination of reference star 

Degrees 


D-22 


PHILCO 



SPACE & RE-ENTRY SYSTEMS DIVISION 
Philco-Ford Corporation 



TR-DA2154 


Table 10 Equation of Motion Parameters 


Index 

Column 

Quantity 

Units 


1 

4 means "equation of motion parameter" 



2-7 

Mnemonic message 


1 

8-9 

Astronomical Unit 


2 

10-11 

Earth's Mass 

Fraction of Sun's Mass 

3 

12-13 

Moon’s Mass 

Fraction of Sun's Mass 

4 

14- 15 

Venus ’"Mass 

Fraction of Sun's Mass 

5 

16-17 

Mars 1 Mass 

Fraction of Sun's Mass 

6 

18-19 

Jupiter's Mass 

Fraction of Sun's Mass 

7 

20-21 

Saturn's Mass 

Fraction of Sun's Mass 

8 

22-23 

Mercury's Mass (not included) 

Fraction of Sun's Mass 

9 

24-25 

Second Zonal Harmonic 

Dimensionless 

10 

26-27 

Third Zonal Harmonic 

Dimensionless 

11 

28-29 

Fourth Zonal Harmonic 

Dimensionless 

12 

30-31 

Fifth Zonal Harmonic 

Dimensionless 

13 

32-33 

First Longitudinal Harmonic 

Radians 

14 

34-35 

Second Longitudinal Harmonic 

Radians 

15 

36-37 

Third Longitudinal Harmonic 

Radians 

16 

38-39 

Fourth Longitudinal Harmonic 

Radians 

17 

40-41 

First Lunar Gravity Parameter 

(Kilometers)^ 

18 

42-43 

Second Lunar Gravity Parameter 

(Kilometers)^ 

19 

44-45 

Third Lunar Gravity Parameter 

2 

(Kilometers) 

20 

46-47 

First Drag Parameter 

(Kilometers) ^ 

21 

48-49 

Second Drag Parameter 

(Kilometers) * 

22 

50-51 

Solar Radiation Pressure Factor 

•J 

(Kilometers) /(seconds) 

23 

52-53 

Venting Thrust Magnitude 

2 

Kilometers/ (seconds) 

24 

54-55 

Speed of Light Error Factor 

Dimensionless 
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/INPCOM/C- ARRAY VALUES 


LOCATION 



NAME 

COMPILED 

VALUE 

DEFINITION 

i 

Constants | 



i 


1 

KPI 

1.57079632679497 

t 

tt/2 

2 

PI 

3.14159265358979 

rr 

3 

TPI 

j 

6.28318530717959 

2tt 

4 

RTD 

57.2957795130823 

Conversion factor, radians to 
degrees , 

5 

i 

DTR 

.017453292519943 

Conversion factor, degrees to 
radians . 

6 

SPMSD 

86400. 


Seconds per mean solar day. 

7 

RSPMSD 

1.15 740 740740 74E- 5 

Reciprocal seconds per mean 
solar day - 1/86400. 


8 

ASTU 

. 149599E9 

Kilometers per astronomical 




unit. 

9 


299774. 

Speed of light, km/sec. 


10 


Body Constants 


j .10 


11-90 


91-98 


B0DC(I,J) 


SBEV1 constant for step-size formul 


th. 

J body constant for body # I, 
where J = 1 to 8 and I = 1 to 7 
is loaded, per Table 1: B0DC(I,J). 
I = 8 to 10 is available for 
three extra bodies with constants 
J ordered as per Table 1 (page 
10.) 


[Initial Conditions! 
SECO 

TARG 


0.0 


5.0 


Starting time, seconds from 
epoch. 

Target body number, (Mars). 


* Unused cells. 
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r 


LOCATION 

NAME 

COMPILED 

VALUE 

DEFINITION' 

113 

PI(2) 

100. 

Diagonal elements P(l,l) 

114 

PI(3) 

100. 

P<2,2) 

115 

PI(4) 

100. 

P(3 ,3) 

116 

PI(5) 

.0001 

P(4,4) 

117 

PI(6) 

.0001 

P(5,5) 

118 

PI (-7) 

.0001 

P(6,6) 

119-123 

• 

PI(8- 12) 

0.0 

From left to right starting 1 
beyond the diagonal: the remain- 

ing 5 elements of row 1. 

124-127 

PIC13-16) 

0.0 

The remaining 4 elements of row 
2. 

128-130 

PI(17-19) 

0.0 

Remaining 3 elements of row 3. 

131-132 

PI(20-21) 

0.0 

Remaining 2 elements of row 4. 

133 

PI(22) 

0.0 

Remaining element of row 5. 

Guidance Spe 

cif ication 



134- 155 

PARI 


Specifications for the symmetric 
initial PAR covariance matrix. 

134 

PARI(l) 

-1. 

Type of coordinates to be des- 
cribed in PARI (2 thru 22), -1. 
to signify the initial PAR is 
identical to the initial upper 



- 

left 6x6 of P. (The using program 
would thus ignore PARI (2 thru 22)). 

156 

PRED 

0.0 

Prediction key (no prediction for 
zero or negative value unless 
guidance is included). 

157 

GUID 

0.0 

Guidance law; set 
< 0, for no guidance 

1. fixed time of arrival 

2. constant target energy 
> 3. minimum energy. 
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f 

LOCATION j 

NAME 

COMPILED 

VALUE 

' 1 1 1 

i 

DEFINITION • ; 

158-162 j 

1 

GUIDT 


Five chronologically ordered 
guidance times referenced to 
starting epoch, in format: days x 
10^ + hours + minutes x 10“ 2 + 
seconds x 10 - ^. 

158 

! 

GUIDT(l) 

1000. 

' 

10 days from epoch or j 

Feb 20, 1975 at 2 h 1® 25.146 s . j 

159 

GUIDT (2) 

10000. 

100 days from epoch. 

160 

GUIDT (3) 

15000. 

150 days. 

161 

j 

GUIDT (4) 

20000. 

200 days. 

162 

GUIDT (5) 

25000. 

250 days, (never to be executed 
for the above Earth -£lars sample 
because actual flight time is 
around 235 days) . 

163 

(EXER) 

30000. 

Percentage error for monitoring . 
guidance correction. The 
numerical value loaded by this 
block data is erroneous. The 
user should overlay an appro- 
priate value when using the 
guidance option. 

164 

GUIDI(l) 

10. 

Resolution error standard devia- 
tion, (10 meters/second) . 

165 

GUIDIC2) 

1. 

Proportional error standard 
deviation, (1 %) . 

166 

GUIDI (3 ) 

1. 

Pointing error standard devi- 
ation, (1 degree). 

167-584 


0.0 

Certain of the following arrays 
are herein set to zero for initial- 
ization. 

167-196 

XTRB 

0.0 

XTRB(I,J) to contain 10 data values^ 
I = 1 to 10, for up to three extra 
bodies ’J = 1 to 3'. H is date, f date, 
body center, coordinate type, and 6 
coordinates, same in ‘order and units 
as C(102 thru .111) . (XTRB(3 ,J)=0. ; 
signals no extra bodies included.) 
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} LOCATION 
\ 

NAME 

COMPILED 

VALUE 

DEFINITION 

jl97-199* 



• 

200-475 

S 


S(I,J) to contain 23 data values, 
I = 1 to 23 for each station # J, 
J = 1 to 12. 

476-566 



B(I) to contain up to 91 beacon 
data values. 

567-584 

OB 


OB (I) to contain up to 18 onboard 
data values 

585-608 

EMP 

1. 

The EMP array serves 3 purposes in 
the Mark II Error Propagation Pro- 
gram, whereas in the Patched Conic 
Program only the first purpose is 
nominally served. 




1. An ordered array of variances 
on equation of motion error sources 
for use in error propagation. 

(The Patched Conic Program is un- 
able to execute equation of motion 
error propagation.) 




2. An ordered array of keys for 
subroutine VAREQ, where zero values 
set logic for Subroutine FBUS to 




make certain omissions; e.g. , 
omissions of related perturbations 
in the gradient calculation. (ERP 
portion of Mark II) . 




3. An ordered array of keys for 
subroutine MPSENS, where for each 



• 

EMP (I) = 0.0 the sensitivities of 
the state with respect to equation 
of motion error source I are com- 
puted. (PINT portion of Mark II). 

609 

ZH(1) 

- . 1082E-2 

| ZHi are zonal harmonic coefficients, 

610 

ZH(2) 

-.23E-5 

1 ^i+l’o Eart h’ s gravitational 

611 

ZH(3) 

-.18E-5 

j field (dimensionless). 

612 

ZH(4) 

0 . 

J 
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LOCATION 

NAME 

COMPILED 

VALUE 

DEFINITION 

513-622 

TH 


Tesseral harmonics where: 

J is the m c h longitudinal 

n,ia th 

harmonic coefficient of the n — 

order, (no units); and X is 

n,m 

the geographic longitude reference 

of J , where -rr < X < tt, 

n,m nm — 

radians . 

613 

TH(1) 

0. 

J 21 

614 

TH(2) 

0. 

*21 

615 

■ TH(3) 

-.120E-5 

J 22 

616 

TH(4) 

-.460 

*22 

617 

TH(5) 

-1.9E-6 

J 31 

618 

TH(6) 

• .080 

*31 

619 

TH(7) 

- . 14E-6 

J 32 

620 

TH(8) 

-.293 

*32 

621 

TH(9) 

- . 10E- 6 

J 33 

622 

TH(10) 

.743 

*33 

623 

ABC 


, Mass- normalized moments of inertia 
of the moon used for computing 
the triaxial gravity perturbation, 
(km 2 ) . 

623 

ABC(l) 

1.20926K6 

a 

624 

ABC (2) 

1.20952E6 

b 

625 

ABC (3) 

1.21003E6 

c 

626 

DRAGC(l) 

. 109 • 

Atmospheric drag coefficients of 
the model: 

> 627 

DRAGC(2) 

. 14294 

-C^e ^2^1 V 1 V, (f>oth have units 
of Km *.) 

628 

TMAG 

0. 

Magnitude of^thrusting accelera- 
tion (km/sec ) . 

629 

SPK 

.94e7 

Solar radiation pressure coeffi- 
cient (km /sec 2 ). 
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LOCATION NAME 


COMPILED 

VALUE 


DEFINITION 


DUTD 


.4050926E-3 


Discrepancy between Ephemeris Time 
and Universal Time, (days). 

The remainder of INPCOM are herein 
set to zero for initialization 
of the following data. 


ICASE 


HEAD(1-12) I 0.0 


HEAD(13-24) 


Case counter for trajectories 
generated and written on binary 
tape by the PINT 'portion of the 
Mark II Program. (ICASE is 
logically identical to KASE, see 
C(674) . ** 

Case counter for ERP cases. 

Case counter for trajectories 
generated and written on binary 
tape by Subroutine CONW. Each 
trajectory generated is numbered 
consecutively and written on tape, 
so that subsequent tape-using 
routines may select the trajectory 
desired by examining the tape for 
this identification number. 
Initially , the case counter must 
be zero to pr op er ly position the- 
binary tape. 

12 word alphanumeric header on 
trajectory tape, # 10. -This 
header is written by CONW or PINT, 
and used by ERP to identify the 
trajectory tape required. HEAD 
is required input. 

12 word alphanumeric header 

1. Describing the ERP case being 
computed . 

2. Written on the special output 
tape if ITAPE =0.0. 

3. Used by SPOUT to identify the 
special output tape required. 


** Refer to footnote on following page. 
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ABM 


COMPILED 

VALUE 

DEFINITION 

699 

ITAPE 

0.0 

ITAPE “ 0, causes special output 
tape # 12 to be written during 
execution of ERP. 

700 

IQCAS 

0 

Case counter for cases written 
on special output tape. Sub- 
sequently read by SPOUT to identify 
the desired case to process,* 


** Program' logic requires that ICASE, XKAS, KASE, and IOGAS be initially 
2 ero, as provided herein. These case counters are automatically 
incremented by the Mark II Program and should never be altered by 
an alternate BLOCK DATA or by ROVLEY .input . 
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TABLE 1: BODC(I , J) 


1 Name 


EARTH 


VENUS MARS 


SATURN 


JUPITR 


gravity 
2 constants, 

, 3 / 2 

kin /sec 


398603.2 


4900.7588 


•13271545E12 324769.5 42977.8 


3.791870E7 1.267106E8 


RPL(I) , 

3 semi-major 
axis, km 


6378. 165 


1738. 


695500. 


6200. 3400. 


60400. 


71350. 


Semi-minor 
axis , km 


6356.5838 4 1738. 


695500. 


54050. 


66600 


RPAT(I) 
sphere of 
influence , 
km 


925000. 


66000. 


1.E10 


616000. 565000. 


. 546E8 


Rotation 
6 rate .rad/sec 


. 72921152E-4 . 266169952E- 5 0.0 


, 70882177E-4 .170553347E-3 . 177491110E-3 


Max step 
7 size, 
seconds 


43200. 


21600. 


864000. 


8 Interpolation 40. 
interval, days 


43200. 43200. 


86400. 


172800. 
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D.2.3 SPOUT 


Input for a SPOUT case consists of: 

1. EXEC card (3 in column 5). 

2. Header for special output tape and IC data, ending with the 
blank. (Header is type K = 0) , 

Blank card must appear. 
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Input Data For the SPOUT Option 


Nominal Value Program Name Location 

Description 

1 


IC(1) 

Case (from tape) to use for output. 

1 


IC(2) 

Regular output point number (counter 
IKY) at which to start outputting. 

241 

KEND 

IC(3) 

Regular output point number at which 
to end output. 

-10 

KDEL 

IC(4) 

Output interval (SPOUT will output 
every KDEL points of type KTYP from 
IC(2) to KEND). 

2 

IPLT 

IC (5) 

Plotting interval (relative to 
IC(4) ) . 

5 

NTYPE 

IC(6) 

No. of output types (This should not 
be changed unless more coding is 
added) . 

1 

KTYP 

IC(7) 

Critical event type which will control 
output . 


ITBLE 


IC(7) ** 0, time control 

1, regular output point 
control 

K, type K event control (See 
KEY in description of tape) 

This is the output choice table. 

To choose a type of output, set: 

1 


IC(10) 

= 1 to get time and state. Details 
of output are set in IC(20). 

1 


IC(ll) 

= 1 to get matrix output. Details of 
output are set in IC(40) . 

0 


IC(12) 

(Space reserved for measurement out- 
puts . ) 

1 


IC(13) 

= 1 to plot. Selection of variables 
to be plotted is in IC(60). 

0 


IC(14) 

(Space reserved for guidance outputs.) 

1, 2, 

3, 0 

10(20)'! 

Details of state output to be printed. 

0, 1, 

111, 0 

IC(30) J 

Up to ten different combinations of 
body center and coordinate system may 
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Nominal Value Program Name Location 


IC(40) 


Description 


be chosen. The vehicle position and 
velocity are output in the coordinate 
system specified in the IC(30) list 
centered at the body set in "the 
corresponding slot of the IC(20) 
list. A zero in the IC(20) list 
terminates the state outputs. The 
body center numbers are; 


1 Earth 

2 Moon 

3 Sun 

4 Venus 

5 Mars 

6 Saturn 

7 Jupiter 


The coordinate system code is made up 
by summing: 


0 

(1950) 


or 1 

(date) 


•f 


0 10 
(equator) (ecliptic) 


or 20 

(selenographic) 


+ 


0 100 or 200 • 

(Cartesian) (spherical) (elements) 


Matrix output options. Zero terminates 
the list. Each entry is made up 
according to the 'code: 

0 10 or 20 

(inertial) (N.V.W. ) (elements) 

4 * 

1 or 2 

(matrix + STD) (Eigen values + EUler) 
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Nominal 

0 

0 

0 

0 

0 

0 

1 

0 


25 

-26 

-25 

-26 

107 

0 


0.0 

0.0 

0.0 

0.0 


Value Program Name Location 
ITEMP IC(50) 


1C (60) 


TSTART 

WC(1) 

TEND 

WC(2) 1 

DELT 

WC(3) j 

TT0L 

WC(4) ^ 


Description 

This is an array which allows 
output at any critical event 
of a given type, even though 
the run is regularly searching 
for some other±ype, or time. 

Set ITEMP (K)^P0 to get out- 
put at every critical event 
type K, K t IC(7). K is as 
described under record 1, KEY. 
(This is independent of IC(7) 
except that if IC(7) > 1, search 
for events of type K will not 
start until after a record with 
sequence number IC(2). If IC(7) 

£ 1, output at type K events 
will start at the beginning of 
the tape.) 

Plot table. The numbers (in- 
dices) of variables to be plot- 
ted. A zero terminates the 
list. Entries = 25 (RMSP), 

26 (KMSV) , or the MPO number 
of any RMS on the tape (See 
Section 5.1 for description of 
MPO's). Each entry is + for 
normal variable, - for last 
variable of a plot. (For in- 
stance, 25, -26, 25, 0 is the 
same as 25, -26, -25, 0 and 
will give two plots - the first 
of RMSP + RMSV, the second of 
RMSP alone, 

T start 

only appli- T end 
cable if , 

IC(7) - 0 * * 

T tolerance (pro- 
gram will accept 
any value within 
At of T^ as the 
value being searched 
for). 
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D.2.4 CONW 

Input for a CCNW case consists of: 

1) EXEC card (4 in column 5). 

2) C-array data and headers of both types, ending with 
a blank.- 

The b h v>.k card must appear. 


Location 

C<99) 

C(100) 

C(101) 


C(102) 

C(103) 

C(104) 

C(105) 


CONW Input Data 


Definition 

Zero reference time, from epoch, in the form: 

2 -2 -4 
Days • 10 + hours + minutes *10 + seconds • 10 . 

Target body number. 

Stop time, from epoch, in the same form as zero time. 


Vehicle Initial Conditions 


2 -2 
Year (of 20th century) • 10 + month + day • 10 . 

2 -2 
Hours « 10 + minutes + seconds • 10 . 

Initial body center number. 

2 

Input type in the form A * 10 + B • 10 + C, where 
A, B and C are interpreted as follows: 

A = 0, Cartesian; B = 0, Equator; C = 0, 1950 Epoch; 

A ° 1, Spherical; B = 1, Ecliptic; C D 1, Epoch of 

cl&tlQ * 

A = 3, Orbital elements; B = 2, Body fixed. ' f 
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C(106-lll) Specify initial state consistent with input type. 


All lengths are in km, all times are in seconds, all 
angles are in degrees. 



Cartes ian 

Spherical 

Orbital Elements 

C(106) 

X 

R 

Semi-major axis 

C(107) 

Y 

Lat 

Eccentricity 

C(108) 

Z 

Long 

True Anomaly 

C(109) 

X 

V 

Long, of ascending node 

C(110) 

• 

Y 

Path 

Angle 

Inclination to X-Y plane 

C(lll) 

Z 

Azimuth 

Arg of periapsis from 
ascending node 


Extra Body Initial Conditions are specified in exactly the same form as 
the vehicle initial conditions.. The input locations are indicated below. 


Location 


Definition 


C(167-176) 

C(177-186) 

C(187-196) 


Extra body one. 
Extra body two. 
Extra body three,. 
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MARK IV : PINT 

Precision Integration 


Subroutine 

Name 

ADOT 

AFTER 

ANTR1 


BACK 

BEGIN 

BLOCK 

DATA 

BLOCK 

DATA 

BLOCK 

DATA 

BODCON 

BUFFIL 

BVEC 

convx 

CROSS 

DATOUT 


Description 

Computes the angle between two given vectors. 

Performs the case calculations which follow trajectory 
integration. 

Provides cartesian components of position and velocity of 
sun, moon, and planets on given date and time. 

Two versions of this subroutine are available: 

1. PANTRY*** computes from mean orbital elements, and 

2. DEPHEM** reads a JPL Ephemeris Tape and interpolates. 

Backspaces a binary tape N logical records or a BCD tape 
N physical records . 

Computes cartesian injection state as a function of the 
controls in subroutine REFINE. 

/DQSCON/ Loads data into DQSCON common block for use by 
subroutine DEQS. 

/INPCOM/ Loads data into INPCOM common block for use in the 
Mark IV Program. 

/WCOM/ Loads data into WCOM common block. 


Supplies values of target constraints to subroutine REFINE. 

Locates the required epoch on the JPL Ephemeris Tape and 
reads the appropriate data for use in subroutine DEPHEM. 

Computes the miss-vector components of the orbit, relative 
to the target body, from the cartesian state. 

Converts input state to cartesian, equatorial, 1950 system 
and outputs results. 

Computes the vector cross product." 

Computes and outputs calendar and Julian dates. 


* Refer to the bracketed subroutine for further description; e.g., DOT 
is described in the ADOT subroutine writeup. 

** Not required by approximate ephemeris package, ANTRl- PANTRY. 

*** Not required by JPL Tape Ephemeris package, ANTR1-DEPHEM. 
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MARK IV: PINT 


Subroutine 

Name 

DEPHEM** 


DEQS 

DERIV 

DOT 

DRAG 

EHA 

EL2EX*** 

ENCKE 

EQTOR 

ERROUT 

EXCOV 

EXIN 

EXPAND*** 

FIEF 

FIFL 

FIST 

FNDMXN 

FNORM 


Description 

Computes position and velocity of planets and moon at any 
Julian date, using Everett's formula to interpolate from a 
JPL Ephemeris Tape. 

Integrates n first or second order differential equations. 

Computes the second derivatives of the sensitivities of the 
state to motion parameters. 

Computes the vector dot product, (ADOT*) . 

Computes acceleration due to atmospheric drag, the gradient 
thereof, and partial derivatives of same wrt the drag 
coefficients . 

Computes the Greenwich hour angle and the earth's angular 
ve locity. 

Converts mean orbital elements to cartesian position and 
velocity. 

Calculates the inverse- square perturbing acceleration due 
to a central body. 

Computes the transformation from mean equator and equinox of 
1950.0 to mean equator and equinox of date; computes mean 
obliquity of date. 

Provides programmed response to anticipated errors. 

Directs the flow of the entire Mark IV Error Propagation 
Program. 

Reads interpolation coefficients from binary tape and 
reconstructs an N dimensional vector X at time T. 

Solves Kepler's equation by series in terms of eccentricity 
and mean anomaly. 

Spaces a binary tape forward over N files, or backward over 
N + 1 files . 

Locates a specified case on binary tape, and reads Record 1 
of this case. 

Rewinds the binary input tape, and reads and checks the 
header. 

Finds the maximum or minimum of a function. 

Computes the magnitude of a vector, (ADOT*). 
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MARK IV: PINT 


Subroutine 

Name 

FSUB 

GOTOR 

GRAVTY 

GTRAN 

GTRN 

GTSN 

INDEX 

INTCOF** 
' INVERT 
MISS1 

MOON 

MPSENS 

MTMPLY 

MVTRN 

NUTATE 


ORB 

0RB2X 

ORIENT 

OSUB 


Description 

Computes the second derivatives of the perturbing accelera- 
tions, the integrated transition matrix, and stopping 
functions for trajectory integration. 

Solves Kepler's equation. 

Computes the acceleration due to the central body's 
gravitational field. 

Computes the coordinate transformation from equator and 
equinox of 1950 to selenographic or true equator and 
Greenwich of date. 

Generates one or a sequence of rotational transformations 
from input angles. 

Generates one or a sequence of rotational transformations 
from input sines and cosines. 

Sets initial controls, constraint tolerances, and indicators 
for refining the trajectory. 

Computes interpolation coefficients for subroutine DEFHEM. 
Inverts a matrix. 

Computes and outputs the matrix of partials at the target 
vector wrt the state at time t, (normally the end point). 

Computes the perturbing accelerations (and gradient thereof) 
caused by the triaxiality of the lunar gravity. 

Initializes integration of Motion Parameter Sensitivities 
from a vehicle ephemeris tape. 

Forms the product of any two matrices . 

Computes the product of a 3 x 3 matrix and a 3 x N matrix. 

Computes the transformations from earth's mean equator, 
equinox to earth's true equator, equinox and/or moon's true 
equator, node. 

Computes and outputs orbital elements. 

Converts orbital elements to cartesian position and velocity. 

Finds the angle through which one vector must be rotated about 
a second vector in order that the first vector should form 
a given angle with a third vector. 

Prints trajectory information, writes a binary tape of the 
trajectory, and calls for a center shift at patch points. 
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MARK IV: PINT 


Subroutine 

Name 

OUTX 

OWRITE 

PATCH 

PERT 

PINT 

PSTART 

REFINE 

REST 

ROTAIT 

ROVLEY 

rvan’ 

RVOUT 

SEVENL 

SHIFT 1 

SIXA 

SIXB 

SIXL 


SOLARP 

SPER 


Description 

Writes out cartesian and spherical position and velocity 
components . 

Writes output and spaces tape during integration for motion 
parameter sensitivities. 

Calculates a constraint function for lunar and planetary 
approach. 

Computes the n-body perturbing acceleration and gradient 
thereof . 

Controls the generation of precision integration trajectories. 

Initializes' the Mark IV Error Propagation Program for the 
precision integration option. 

Iterates on a set of launch controls until the resulting 
trajectory satisfies a prescribed set of target arrival 
conditions . 

Performs the trajectory shifting required at rectification. 
Rotates two vectors in a plane. 

Reads fixed, floating, and alphanumeric input data into core. 

Converts spherical coordinates to cartesian. 

Converts cartesian position and velocity to spherical 
coordinates . 

Calls subroutine PSTART or AFTER. 

(SEVENL is the PROGRAM for the (7,0) primary overlay.) 

Calculates position and velocity relative to ephemeris bodies 
and extra bodies. 

Calls subroutine TRAJ. 

(SIXA is the PROGRAM for the (6,1) secondary overlay.) 

Calls subroutine REFINE. 

(SIXB is the PROGRAM for the (6,2) secondary overlay.) 

Calls one of 3 secondary overlays. The (6,0) primary overlay 
makes the DEQS precision integration package available to 
TRAJ or REFINE, (or to PLNTDQ for SEARCH). 

(SIXL is the PROGRAM for the (6,0) primary overlay.) 

Computes the acceleration due to solar radiation pressure and 
the gradient thereof, or the partial derivative of same. 

Computes spherical coordinates of a cartesian 3-vector. 
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Subroutine 

Name 

SSIZE 

■ STEPD 

STEPI 

STEP! 

TARGT1 

TCONIC 

TFRAC 

THRUST 

TIMEC 

TIMED 

TIMES 

TPFL 

TPST 

TRAJ 

UPDATE*** 

VNORM 

XOUT 


MARK IV: PINT 


Description 

Computes a starting integration step size. 

Computes new cartesian state on a conic, given the old state 
plus incremental time or true anomaly. 

Computes the array of orbital elements to define a conic for 
subroutine STEPT. 

Computes cartesian state at given time on given conic. 

Computes the transformation matrix to target coordinates. 

Calculates conic time as a function of true anomaly. 

Updates time in whole and fractional days from epoch. 

Calculates the acceleration due to thrust, the gradient 
thereof, or the partial derivative of same. 

Converts calendar date and time to whole and fractional 
days from January 1, 1950. 

Converts time from (days, hours, minutes, seconds) to 
seconds . 

Converts time from seconds to alphanumeric days, hours, min- 
utes , and seconds . 

Writes Record 1 of each case on binary tape. 

Rewinds a binary tape, writes a header, then an End of File 
on this tape. 

Drives the integration for 1) regular trajectories or tape- 
write, 2) computing equation of motion sensitivities. 

Updates mean orbital elements in time from one epoch to 
another . 

Normalizes a vector and also computes the magnitude, (ADOT*) . 

Computes interpolation coefficients and writes them on a 
binary tape. 
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MARK XV: ERP 

Error Propagation 


Subroutine 

Name 


Description 


ANTR1 


BACK 

BCHNG 

BEACH 

BLDP1 

BLOCK 

DATA 

BLOCK 

DATA 

BUFF XL** 
CHNG 
CONVP 
CRITA 


Provides cartesian components of position and velocity of, 
sun, moon, and planets on given date and time. 

Two versions of this subroutine are available: 

1. PANTRY*** computes from mean orbital elements, and 

2. DEPHEM** reades a JPL Ephemeris Tape and interpolates. 

Back spaces a binary tape N logical records or a BCD tape 
N physical records. 

Makes beacon measurements and appropriately updates the 
expanded covariance matrix. 

Supplies information pertaining to the next beacon critical 
event . 

Builds the expanded P covariance matrix and outputs the error 
sources included. 

/INPCOM/ Loads data into INPCOM common block for use in the 
Mark IV Program. 

/DQSCON/ Loads data into DQSCON common block for use by sub- 
routine DEQS. 

Locates the required epoch on the JPL Ephemeris tape and 
reads the appropriate data for use in Subroutine DEPHEM. 

Makes earth-based tracking station measurements and appropri- 
ately updates the expanded covariance matrix. 

Outputs the input covariance matrix and converts it to a 
6x6 in equator of 1950.0. 

Outputs when a body starts and stops occulting the vehicle. 


CRITO 

CROSS 


Outputs station and beacon in-view, out-of-view critical 
events . 

Computes the vector cross product. 


* Refer to the bracketed subroutine for further description; e.g., DOT 
is described in the ADOT subroutine writeup. 

** Not required by approximate ephemeris package, ANTR1-PANTRY. 

*** Not required by JPL Tape Ephemeris package, ANTR1-DEPHEM. 
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MARK IV: ERP 

Error Propagation 


Subroutine 

Name 


Description 


DATOUT 

DEPHEM** 

DEQS 

DOT 

DRAG 

EHA 

EIGEN 

ELLINT 

EL2EX*** 

EQTOR 

ERP 

ERROUT 

EV.DEL 

EXCOV 

EXIN 

EXPAND*** 

FBUS 

FIEF 

FIFL 


Computes and outputs calendar and Julian dates. 

Computes position and velocity of planets and moon at any 
Julian date, using Everett’s formula to interpolate from a 
JPL Ephemeris Tape. 

Integrates n first or second order differential equations. 
Computes the vector dot product, (ADOT*). 

Computes acceleration due to atmospheric drag, the gradient 
thereof, and partial derivatives of same wrt the drag 
coefficients. 

Computes the Greenwich hour angle and earth's rate. 

Computes eigenvalues and eigenvectors of a real symmetric 
matrix. 

Evaluates the incomplete elliptic integrals of the first 
and second kinds . 

Converts mean orbital elements to cartesian position and 
velocity. 

Computes the transformation from mean equator and equinox of 
1950.0 to mean equator and equinox of date; computes mean 
obliquity of' date. 

Controls the logical flow of error propagation computations. 

Provides programmed response to anticipated errors. 

Orders the critical events chronologically and determines 
the advancing step size. 

Directs the flow of the entire Mark IV Error Propagation 
Program. 

Reads interpolation coefficients from binary tape and 
reconstructs an N dimensional vector X at time t. 

Solves Kepler's equation by series in terms of eccentricity 
and mean anomaly. 

Supplies the DEQS integration routine with the derivatives 
to integrate the variational equations. 

Spaces a binary tape forward over N files, or backward over 
N + 1 files. 

Locates a specified case on a binary tape, and reads record 
1 of this case. 
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Subroutine 

Name 

FIST 

FNORM 

FOURA 

FOURB 

FOURG 

FOURD 

FOURE 

FOURL 


GETEAN 

gmiss 

GORE 

GOTOR 

GRAVTY 

GTRN 

GTSN 

INTCOF** 

INVERT 

IAYOVR 


Description 

Rewinds the binary input tape , and reads and checks the 
header. 

Computes the magnitude of a vector , (ADOT*) . 

Calls subroutine ONBRD at onboard measurement critical 
event, or calls QiESS at body center change. 

(FOURA is the PROGRAM for the (4,1) secondary overlay.) 

Calls subroutine GORE at a guidance correction critical 
event . 

(FOURB is the PROGRAM for the (4,2) secondary overlay.) 

Calls subroutine CHNG as required at a station measurement 
critical event. 

(FOURC is the PROGRAM for the (4,3) secondary overlay.) 

Calls subroutine ONBRAD or BCHNG at a critical event requir- 
ing onboard radar or beacon measurements ,. respectively. 
(FOURD is the PROGRAM for the (4,4) secondary overlay.) 

Calls subroutine DEQS when required by subroutine PHIX. 
(FOURE is the PROGRAM for the (4,5) secondary overlay.) 

Performs error propagation and calls subroutine 0UT1 to 
output results at critical events. 

(FOURL is the PROGRAM for the (4,0) primary overlay). 

Computes the transformation from EE50 coordinates to body- 
fixed coordinates of instant. 

Computes the guidance sensitivity matrix from a specified 
patch point to the end point. , 

Makes a guidance correction. 

Solves Kepler's equation. 

Computes the acceleration due to the central body's gravi- 
tational field. 

Generates one or a sequence of rotational transformations 
from input angles. 

Generates one or a sequence of rotational transformations 
from input sines and cosines. 

Computes interpolation coefficients for subroutine DEPHEM. 
Inverts a matrix. 

Reads station, beacon, onboard, and equation of motion data 
cards and converts to proper units . 
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MARK IV: ERP 

Error Propagation 


Subroutine 

Name 

MOON 

MTMPLY 

MTRX 

MVTRN 

NUTATE 

OBUS 

ONBRAD 

ONBRD ' 

ONEL 

ORB 

OUM 

OUTP 

OUT1 

PARAB 

PCHNG 

PERT 

PHIX 

PHIZ 

PRESTO 


Description 

Computes the perturbing accelerations (and gradient thereof) 
caused by the triaxiality of the lunar gravity. 

Forms the product of any two matrices. 

Multiplies a 6 x 6 matrix times the upper left 6 x 6 of a 
given matrix of equal or larger dimension. 

Computes the product of a 3 x 3 matrix and a 3 x N matrix. 

Computes the transformation(s) from earth's mean equator and 
equinox to earth's true equator and equinox, and/or to moon's 
true equator and node. 

Prints a message when called by DEQS. 

Makes onboard radar measurements and appropriately updates 
the expanded covariance matrix. 

Makes onboard optical measurements and appropriately updates 
the expanded covariance matrix. 

Calls subroutine PRESTO and, if necessary, GMISS to initialize 
the trajectory and other parameters for error propagation. 
(ONEL is the PROGRAM for the (1,0) primary overlay.) 

Computes and outputs orbital elements. 

Outputs the RMS uncertainty in the miss vector. 

Outputs the P and Par covariance matrices in Darboux coordi- 
nates . 

Outputs the state, calculates and outputs RMS values and 
target miss vector on output tape and binary tape. 

Fits a parabola through three points. 

Updates the expanded covariance matrix at an observation. 

Computes the n-body perturbing acceleration and gradient 
thereof . 

Computes the state transition matrix from state at current 
time to state at specified time, and calls FBUS if equation 
of motion error sources are included. 

Computes the state transition matrix (on one conic). 

Reads overlay input for common and initializes the tape 
read and various parameters for the Error Propagation. 
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MARK IV: ERP 

Error Propagation 


Subroutine 

Name 

PROP 

PUTIN 

QUARTO 


Description 

Controls the updating of the state vector and the covariance 
matrix. 

Reads in measurement error source keys. 

Finds the solutions to the quadratic equation. 


ROVLEY 

ROYAL 

RTIMS 

SBEV1 

SCOT 

SHIF2 

SHUFLP 

SOIARP 

SORDR 

SPER 

STABEC 

STASH 

STAT 

STATP 

STEPI 

STEPT 

ICONIC 


Reads fixed, floating, and alphanumeric input data into core. 

Outputs input random, bias, and time errors associated with 
station and beacon measurement error sources. 

Controls input of control times and changes in measurement 
treatment . 

Computes vehicle in-view and out-of-view critical events for 
earth-based tracking stations or beacons. 

Finds the expected value of a given matrix. 

Calculates position and velocity relative to ephemeris bodies 
and extra bodies . 

Rearranges the expanded covariance matrix when error sources 
are added, deleted, or considered differently. 

Computes the. acceleration due to solar radiation pressure and 
the gradient thereof, or the partial derivative of same. 

Sorts an array X in ascending order, while preserving the 
correspondence between array X and array NX. 

Computes spherical coordinates of a cartesian 3-vector. 

Sets up logic and loads critical event arrays for observations 
made by stations, beacons, and on board optical. 

Supplies information for next earth-based tracking station 
critical event. 

Computes inertial coordinates of a tracking station and the 
orthogonal transformation relating inertial cartesian plane 
to local tangent plane North-East-Down. 

Calculates the partials of earth-based tracking station 
measurements with respect to the extended state vector. 

Computes an array of orbital elements to define a conic for 
subroutine STEPT. 

Computes cartesian state at a given time on a given conic. 
Calculates conic time as a function of true anomaly. 
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Subroutine 

Name 

TFRAC 

THREEL 

THRUST 

TIMED 

TIMES 

TPFL 

TPST 

TRANP 

TRDB 

TWOL 


update*** 

UPPT 

VAREQ 

VENT 

VNORM 


MARK IV: ERP 

Error Propagation 


Description 

Updates time in whole and fractional days from epoch. 

Galls subroutine SBEV1 to compute on-off critical events for 
stations or beacons, as required. 

(THREEL is the PROGRAM for the (3,0) primary overlay). 

Calculates the acceleration due to thrust, the gradient 
thereof, or the partial derivative of same. 

Converts time from (days, hours, minutes, seconds) to 
seconds . 

Converts time from seconds to alphanumeric days, hours, 
minutes, and seconds. 

Writes Record 1 of each case on binary tape. 

Rewinds a binary tape, writes a header, then an End of File 
on this tape. 

Transforms the covariance matrix from one inertial frame to 
another via a given transformation matrix. 

Computes Darboux or local tangent plane transformation matrix. 

Calls subroutine RTIMS to read in and organize measurement 
requirements, and/or calls OUTP to print the covariance 
matrix. 

(TWOL is the PROGRAM for the (2,0) primary overlay). 

Updates mean orbital elements in time from one epoch to 
another . 

Propagates the expanded covariance matrix in time. 

Performs initializations for the integration of variational 
equations when equation of motion error sources are 
considered. 

Initializes the KEV and EVNT critical event arrays and other 
parameters . 

Normalizes a vector and also computes the magnitude, (ADOT*) . 


E-ll 


PHILCO 



SPACE £ RE-ENTRY SYSTEMS DIVISION 
Philco-Ford Corporation 



TR-DA2154 


MARK IV : SPOUT 

Special Output 


Subroutine 

Name 

ADOT 

ANTR1 


BLOCK 

DATA 

BLOCK 

DATA 

BUFF I L** 

COVOUT 

CROSS 

DATOUT 

DEPHEM** 

DOT 

EIGEN 

EIGHTL 

EL2EX*** 

EQTOR 


Description 

Computes the angle between two given vectors. 

Provides cartesian components of position and velocity of 
sun, moon, and planets on given date and time. 

Two versions of this subroutine are available: 

1. PANTRY*** computes from mean orbital elements, and 

2. DEPHEM** reads a JPL Ephemeris Tape and interpolates. 

/INPCQM/ Loads data into INPCOM common block for use in the 
MARK IV Program. 

/RVR/ Loads data into RVR common block. 

Locates the required epoch on the JPL Ephemeris Tape and 
reads the appropriate data for use in Subroutine DEPHEM, 

Outputs the covariance matrix as requested by options in the 
special output program. 

Computes the vector cross product. 

Computes and outputs calendar and Julian dates. 

Computes position and velocity of planets and moon at any 
Julian date, using Everett's formula to interpolate from a 
JPL Ephemeris Tape. 

Computes the vector dot product, (ADOT*). 

Computes eigenvalues and eigenvectors of a real symmetric 
matrix- 

Calls subroutine SPOUT. 

(EIGHTL is the PROGRAM for the (8,0) primary overlay.) 

Converts mean orbital elements to cartesian position and 
velocity. 

Computes the transformation from mean equator and equinox 
of 1950.0 to mean equator and equinox of date; computes 
mean obliquity of date. 


* Refer to the bracketed subroutine for further description; e.g., DOT 
is described in the ADOT subroutine writeup. 

** Not required by approximate ephemeris package, ANTE 1-PANTRY . 

*** Not required by JPL Tape Ephemeris package, ANTRl- DEPHEM. 
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MARK IV: SPOUT 


Subroutine 
Name 


Description 


ERROUT 

EXCOV 

EXPAND*** 

EXTRA 

FIEF 

FIFL 

FINFO 

FNORM 

GTRN 

GTSN 

INTCOF** 

LOUT 

MEAS 

MVTRN 

NMOUT 

NUTATE 


PLOTZ 

PLT1M 


Provides programmed response to anticipated errors. 

Directs the flow of the entire MARK IV Error Propagation 
Program. 

Solves Kepler's equation by series in terms of eccentricity 
and mean anomaly. 

Reads the extended output binary tape produced by an ERP run, 
processes this data according to input options, and saves 
this data for plotting. 

Spaces a binary tape forward over N files , or backward over 
N + 1 files. 

Locates a specified case on the extended output binary tape, 
and reads Record 1 of this case. 

Locates a specified output set on the extended output binary 
tape. 

Computes the magnitude of a vector, (ADOT*) . 

Generates one or a sequence of rotational transformations 
from input angles . 

Generates one or a sequence of rotational transformations 
from input sines and cosines. 

Computes interpolation coefficients for subroutine DEPHEM. 

Normalizes and outputs the P covariance matrix. 

Outputs measurements and associated RMS quantities. 

Computes the product of a 3x3 matrix and a 3xN matrix. 

Outputs the vehicle state in the coordinate system(s) 
requested by input. 

Computes the transformation^ ) from earth's mean equator 
and equinox to earth's true equator and equinox, and/or to 
moon's true equator and node. 

Plots the variables saved during a run of the special 
output program. 

Makes rough plots of selected variables for output by the 
printer. 
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MARK IV : SPOUT 


Subroutine 

Name 

PLTSAV 

RNDLIM 

ROVLEY 

SPER 

SPOUT 

SUBPLT 

TFRAC 

TIMES 

TRANP 

TRDB 


Description 

Accumulates an array of variables for plotting. 

Produces rounded limits for a plot axis. 

Reads fixed, floating, and alphanumeric input data into core. 

Computes spherical coordinates of a cartesian 3-vector. 

Controls the flow of the special output portion of the . 

MARK IV Error Propagation Program. 

Prepares one line of output for the printer plotter sub- 
routine PLT1M. 

Updates time in whole and fractional days from epoch. 

Converts time from seconds to alphanumeric days, hours, 
minutes , and seconds . 

Transforms the covariance matrix from one inertial frame to 
another via a given transformation matrix. 

Computes Darboux or local tangent plane transformation matrix. 


UPDATE*** 

VNORM 

X20RB 


Updates mean orbital elements in time from one epoch to 
another . 

Normalizes a vector and also computes the magnitude, (ADOT*) . 

Computes orbital elements from cartesian positions and 
velocity . 
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MARK IV : CONW 

Patched Conic Tape Generation 


Subroutine 

Name 


Description 


ADOT 

ANTR1 


BLOCK 

DATA 

BUFFIL** 

BVEC 

CONVX 

CONW 

CROSS 

DATOUT 

DEPHEM** 

DOT 

EHA 


Computes the angle between two given vectors. 

Provides cartesian components of position and velocity of 
sun, moon, and planets on given date and time. 

Two versions of this subroutine are available: 

1. PANTRY*** computes from mean orbital elements, and 

2. DEPHEM** reads a JPL Ephemeris Tape and interpolates. 

/INPCOM/ Loads data into INPCOM common block for use in the 
Mark IV Program. 

Locates the required epoch on the JPL Ephemeris Tape and 
reads the appropriate data for use in Subroutine DEPHEM. 

Calculates the target miss vector. 

Converts input state to cartesian, equatorial, 1950 system 
and outputs results . 

Generates a patched conic trajectory and stores it on a 
binary tape for use by other programs. 

Computes the vector cross product. 

Computes and outputs calendar and Julian dates. 

Computes position and velocity of planets and moon at any 
Julian date, using Everett's formula to interpolate from a 
JPL Ephemeris Tape. 

Computes the vector dot product, (ADOT*). 

Computes the Greenwich hour angle and Earth's rate. 


EL2EX*** Converts mean orbital elements to cartesian position and 
velocity. 


* Refer to the bracketed subroutine for further description; e.g., 
DOT is described in the ADOT subroutine writeup. 

** Not required by approximate ephemeris package, ANTRl- PANTRY. 

*** Not required by JPL Tape Ephemeris package, ANTR1-DEPHEM. 
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MARK IV : CONW 


Subroutine 

Name 


Description 


EQTOR 

ERROUT 

EXCOV 

EXPAND*** 

FIVEE 

FNORM 

GOTOR 

GTRN 

GTSN 

INTCOF** 

MISS1 

MVTRN 

NUTATE 

ORB 

0RB2X 

OUTX 

PATCH 

ROTA IT 
ROVLEY 
RVAN 
RVOUT 


Computes the transformation from mean equator and equinox 
of 1950.0 to mean e'quator, equinox of date; computes mean 
obliquity of date. 

Provides programmed response to anticipated errors . 

Directs the flow of the entire Mark IV Error Propagation 
Program. 

Solves Kepler's equation by series in terms of eccentricity 
and mean anomaly. 

Calls CONW to generate a patched conic trajectory. (FIVEL 
is the PROGRAM for the (5,0) primary overlay.) 

Computes the magnitude of a vector, (ADOT*). 

Solves Kepler's equation. 

Generates one or a sequence of rotational transformations 
from input angles. 

Generates one or a sequence of rotational transformations 
from input sines and cosines. 

Computes interpolation coefficients for subroutine DEPHEM. 

Computes and outputs the matrix of partials of_the -target 
vector wrt the state at time t (normally the endpoint). 

Computes the product of a 3x3 matrix and a 3xN matrix. 

Computes the transf ormation(s) from earth's mean equator 
and equinox to earth's true equator and equinox, and/or to 
moon's true equator and node. 

Computes and outputs orbital elements. 

Converts orbital elements to cartesian position and 
velocity. 

Writes out cartesian and spherical position and velocity 
components . 

Calculates a constraint, function for lunar and planetary 
approach. 

Rotates two vectors in a plane. 

Reads fixed, floating, and alphanumeric input data into core. 

Converts spherical coordinates to cartesian. 

Converts cartesian position and velocity to spherical 
coordinates . 
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MARK IV : CONW 


Subroutine 

Name 


Description 


SHIFT 1 

SKET 

SPER 

STEPD 

TARGT1 

TCONIC 

TFRAC 

TIMEC 

TIMED 

TIMES 


Calculates position and velocity relative to ephemeris 
bodies and extra bodies . 

Determines the trajectory type and sets an index accordingly. 

Computes spherical coordinates of a cartesian 3-vector. 

Calculates cartesian state on a conic by stepping time or 
true anomaly. 

Computes the transformation matrix to target coordinates, 
(TARGT*) . 

Calculates conic time as a function of true anomaly. 

Updates time in whole and fractional days from epoch. 

Converts calendar date and time to whole and fractional 
days from January 1, 1950. 

Converts time from (days, hours, minutes, seconds) to seconds. 

Converts time from seconds to alphanumeric days, hours, min- 
utes , and seconds . 


TPST 

TRUEA 

TTGO 

UPDATE*** 

VNORM 

XOUT 


Rewinds a binary tape, writes a header, then an End of File 
on this tape. 

Computes true anomaly from semi-latus rectum, eccentricity, 
and radius in the orbit. 

Computes time-to-go function for patching conics. 

Updates mean orbital elements in time from one epoch to 
another . 

Normalizes a vector and also computes the magnitude , (ADOT*) . 

Computes interpolation coefficients and writes them on a 
binary tape. 
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MARK IV: START-UP 


Subroutine 

Name 


Description 


ADOT 

ANTR1 


BLOCK 

DATA 

BLOCK 

DATA 

BUFFIL** 

CONBR 

CROSS 

DATOUT 

DEPHEM** 


DOT 

EHA 

EL2EX*** 

EQTOR 


Computes the angle between two given vectors . 

Provides cartesian components of position and velocity of 
sun, moon, and planets on given date and time. 

Two versions of this subroutine are available: 

1. PANTRY*** computes from mean orbital elements, and 

2. DEPHEM** reads a JPL Ephemeris Tape and interpolates. 

/INPCOM/ Loads data into INPCOM common block. 

/WCOM/ Loads data into WCOM common block. 

Locates the required epoch on the JPL Ephemeris Tape and 
reads the appropriate data for use in Subroutine DEPHEM. 

Computes the conic trajectory connecting two points in a 
given time. 

Computes the vector cross product. 

Computes and outputs calendar and Julian dates . 

Computes position and velocity of planets and moon at any 
Julian date, using Everett's formula to interpolate from a 
JPL Ephemeris Tape. 

Computes the vector dot product, (ADOT*). 

Computes the Greenwich hour angle and the earth's angular 
velocity. 

Converts mean orbital elements to cartesian position and 
velocity. 

Computes the transformation from mean equator and equinox of 
1950.0 to mean equator and equinox of date. 


* Refer to the bracketed subroutine for further description; e.g., DOT 
is described in the ADOT subroutine writeup. 

** Not required by approximate ephemeris package, ANTR1-PANTRY. 

*** Not required by JPL Tape Ephemeris package, ANTRl-DEPHEM. 
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Subroutine 

Name 

ERROUT 

EXCOV 

EXPAND*** 

FINDV 

FNORM 

GOTOR 

GTRN' 

GTSN- 

INTCOF** 

LAUNCH 

MVTRN 
NINEL . 

NUTATE 

ORB 

ORIENT 

OXJTX 

PLANET 

ROTAIT 

ROVLEY 


MARK IV: START-UP 


Description 

Provides programmed response .to -anticipated errors. 

Directs the flow of the entire MARK IV Error Propagation 
Program. 

Solves Kepler's equation by series in terms of eccentricity 
and mean anomaly. 

Minimizes or finds zeroes of a function of a scalar variable. 
Computes the magnitude of a vector, (ADOT*) . 

Solves Kepler's equation. 

Generates one or a sequence of rotational transformations 
from input angles. 

Generates one or a sequence of rotational transformations 
from input sines and cosines. 

Computes interpolation coefficients for subroutine DEPHEM. 

Computes the control values which join the parking orbit to 
the hyperbolic excess velocity. 

Computes the product of a 3x3 matrix and a 3xN matrix. 

Calls subroutine PLANET then stores solution in WCOM common, 
or calls subroutine SEARCH, 

(NINEL is the PROGRAM for the (9,0) primary overlay). 

Computes the transf ormation(s) from earth's mean equator 
and equinox to earth's true equator and equinox, and/or to 
moon's true equator and node. 

Computes and outputs orbital elements. 

Finds the angle through which one vector must be rotated 
about a second vector in order that the first vector form a 
given angle with a third vector. 

Writes out cartesian and spherical position and velocity 
components . 

Start-up driver - computes approximate interplanetary 
trajectory solutions. 

Rotates two vectors in a plane. 

Reads fixed, floating, and alphanumeric input data into core. 
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MARK IV: START-UP 


Subroutine 

Name 


Description 


RVOUT 

SHIFT 1 

SPER 

START 

STEPD 

TCONIC 

TFRAC 

TIMEC 

TRDB 

TRUEA 

UPDATE*** 

VELASY 

VNORM 


Converts cartesian position and velocity to spherical 
coordinates . 

Calculates position and velocity relative to ephemeris 
bodies and extra bodies . 

Computes spherical coordinates of a 3-vector. 

Computes cartesian injection state as a function of launch 
parameter controls. 

Calculates cartesian state on a conic by stepping time or 
true anomaly. 

Calculates conic time as a function of true anomaly. 

Updates time in whole and fractional days from epoch. 

Converts calendar date and time to whole and fractional 
days from January 1, 1950. 

Computes the transformation to 'local tangent plane or 
Darboux coordinates. 

Computes true anomaly from semi-latus rectum, eccentricity, 
and radius in the orbit. 

Updates mean orbital elements in time from one epoch to 
another. 

Computes the velocity vectors and differences at the end - 
points of the conic section which connects two radii in a 
given time. 

Normalizes a vector and also computes the magnitude, (ADOT*) . 


E-20 


PHILCO 


SPACE S. RE-ENTRY SYSTEMS DIVISION 
Phifeo-Ford Corporation 



TR-DA2154 


MARK IV : SEARCH 


Subroutine 

Name 


Description 


ADOI 

ANTR1 


BLOCK 

DATA 


Computes the angle between two given vectors. 

Provides cartesian components of position and velocity of 
sun, moon, and planets on given date and time. 

Two versions of this subroutine are available: 

1. PANTRY*** computes from mean orbital elements, and 

2. DEPHEM** reads a JPL Ephemeris Tape and interpolates. 

/DQSCON/ Loads data into DQSCON common block for use by sub- 
routine DEQS. 


BLOCK /INPCOM/ Loads data into INPCOM common block. 

DATA 


BLOCK 

DATA 

BUFFIL** 

BVEC 

CONVX 

CROSS 

DATOUT 

DEPHEM** 

DEQS 

DOT 

DRAG 


/WCOM/ Loads data into WCOM common block. 

Locates the required epoch on the JPL Ephemeris Tape and 
reads the appropriate data for use in Subroutine DEPHEM. 

Computes the miss vector components of the orbit, relative to 
the target body, from the cartesian state. 

Converts input state to cartesian, equatorial, 1950 system and 
outputs results. 

Computes the vector cross product. 

Computes and outputs calendar and Julian dates . 

Computes position and velocity of planets and moon at any 
Julian date, using Everett's formula to interpolate from a 
JPL Ephemeris Tape. 

Integrates n first or second order differential equations. 
Computes the vector dot product, (ADOT*). 

Computes acceleration due to atmospheric drag, the gradient 
thereof, and partial derivatives of same wrt the drag 
coefficients . 


* Refer to the bracketed subroutine for further description; e.g., DOT 
is described in the ADOT subroutine writeup. 

** Not required by approximate ephemeris package , ANTR1- PANTRY. 

*** Not required by JPL Tape Ephemeris package, ANTR1-DEPHEM. 
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MARK IV : SEARCH 


Subrout ine 
Name 

EHA 

EL2EX*** 

ENCKE 

ENDCON 

EQTOR 

ERROUT 

EXCOV 

EXPAND*** 

FNDMXN' 

FNORM 

FSUB 

GOTOR 


Description 

Computes the Greenwich hour angle and the earth's rotation 
rate. 

Converts mean orbital elements to cartesian position and 
velocity. 

Calculates the inverse- square perturbing acceleration due 
to a central body. 

Computes the end constraints for SEARCH. 

Computes the transformation from mean equator and equinox 
of 1950.0 to mean equator and equinox of date; computes 
mean obliquity of date. 

Provides programmed response to anticipated errors. 

Directs the flow of the entire Mark IV Error Propagation 
Program. 

Solves Kepler's equation by series in terms of eccentricity 
and mean anomaly. 

Finds the maximum or minimum of a function. 

Computes the magnitude of a vector, (ADOT*) . 

Computes the second derivatives of the perturbing accelera- 
tions, the integrated transition matrix, and stopping 
functions for trajectory integration. 

Solves Kepler's equation. 


GRAVTY 

GTRAN 

GTRN’ 

GTSN 

INTCOF** 

INVERT 

MARFIX 


Computes the acceleration due to the central body's gravita- 
tional field. 

Computes the transformation from equator and equinox of 1950 
to selenographic or true equator and Greenwich of date. 

Generates one or a sequence of rotational transformations 
from input angles. 

Generates one or a sequence of rotational transformations 
from input sines and cosines. 

Computes interpolation coefficients for subroutine DEPHEM. 
Inverts a matrix. 

Computes the transformation from mean equator and equinox of 
1950.0 to mars fixed coordinates. 
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MARK XV: SEARCH 


Subroutine 

Name 

MOON 

MTMPLY 

MVTRN 

NINEL 

NUTATE 

ORB 

ORB2X 

OSUB 

OUTX 

PATCH 

PCON 

PERT 

PHIZ 

PLNTDQ 

PLNTPC 

REST 

ROTAIT 

ROVLEY 

RVAN 

RVOUT 

SEARCH 


Description 

Computes the perturbing accelerations (and gradient thereof) 
caused by the triaxiality of the lunar gravity. 

Forms the product of any two matrices. 

Computes the product of a 3 x 3 matrix and a 3 x N matrix. 

Calls subroutine PLANET then stores results, or calls sub- 
routine SEARCH. 

(NINEL is the PROGRAM for the (9,0) primary overlay.) 

Computes the transformation(s) from earth's mean equator and 
equinox to earth's true equator and equinox, and/or to moon's 
true equator and node. 

Computes and outputs orbital elements. 

Converts orbital elements to cartesian position and velocity. 

Prints integrated trajectory information, writes a binary 
tape of the trajectory. 

Writes out cartesian and spherical position and velocity 
components . 

Calculates a constraint function for lunar and planetary 
approach. 

Computes points on a patched conic trajectory. 

Computes the n-body perturbing acceleration and gradient 
thereof. 

Computes the closed-form conic state transition matrix. 

Drives the integrated trajectory calcuation for SEARCH. 

Drives the patched conic trajectory calculation for SEARCH. 

Performs the trajectory shifting required at rectification 
and patch points. 

Rotates two orthonormal vectors through an angle in their 
mutual plane. 

Reads fixed, floating, and alphanumeric input data into core. 

Converts spherical coordinates' to cartesian. 

Converts cartesian position and velocity to spherical 
coordinates. 

Searches to satisfy selected end constraints by varying 
selected trajectory controls. 
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MARK IV: SEARCH 


Subroutine 

Name 

SETUP 

SHIFT1 

SIXC 

SIXL 


SKET 

SOLARP 

SPER 

SRCHOV 

SSIZE 

START 

STEPD 

STEPI 

STEPT 

TCONIC 

TFRAC 

THRUST 

TIMEC 

TIMED 


Description 

Interprets the input options, sets up the various, indicators 
and parameters required by SEARCH, and prints out requested 
options. 

Calculates position and velocity relative to ephemeris bodies 
and extra bodies. 

Calls subroutine PLNTDQ to provide an integrated trajectory 
for SEARCH. 

(SIXC is the PROGRAM for the (6,3) secondary overlay). 

Calls one of three secondary overlays. The (6,0) primary 
overlay makes the DEQS precision integration package avail- 
able to PLNTDQ for use in SEARCH, or to TRAJ or REFINE as 
required by PINT, 

(SIXL is the PROGRAM for the (6,0) primary overlay.) 

Determines the trajectory type and sets an index accordingly. 

Computes the acceleration due to solar radiation pressure and 
the gradient thereof, or the partial derivative of same. 

Computes spherical coordinates of a cartesian 3-vector. 

Driver for the generalized search capability. 

Computes a starting integration step size. 

Computes cartesian injection state as a function of launch 
parameter controls. 

Computes new cartesian state on a conic, given the old state 
plus incremental time or true anomaly. 

Computes the array of orbital elements to define a .conic for 
subroutine STEPT. 

Computes cartesian state at a given time on a given conic. 

Computes the time from periapsis for a given true anomaly on 
a Kepler ian conic. 

Updates time in whole and fractional days from epoch. 

Calculates the acceleration due to thrust, the gradient thereof, 
or the partial derivative of same. 

Converts calendar date and time to whole and fractional days 
from January 1, 1950. 

Converts a time interval from days, hours, minutes, and 
seconds to seconds. 
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MARK IV : SEARCH 


Subroutine 

Name 

TIMES 


TRDB 

TRUEA 


Description 

Converts a time interval from seconds to days, hours, minutes 
and seconds and sets up an alphanumeric array for subsequent 
output . 

Computes the transformation to Darboux or local tangent plane. 

Computes true anomaly, given semi-latus rectum, eccentricity, 
and radius in the orbit. 


TTGO Computes time-to-go function for patching conics. 

UPDATE*** Updates mean orbital elements in time from one epoch to 
another . 


VNORM Normalizes a vector and also computes the magnitude, (ADOT*). 

XOUT Computes interpolation coefficients and writes them on a 

binary tape. 
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